<!DOCTYPE html PUBLIC "-//W3C//DTD HTML 4.01 Transitional//EN">
<html><head><meta http-equiv="Content-Type" content="text/html; charset=utf-8"><meta http-equiv="X-UA-Compatible" content="IE=edge,IE=9,chrome=1"><meta name="generator" content="MATLAB 2021a"><title>Flux Balance Analysis (FBA) and its variants</title><style type="text/css">.rtcContent { padding: 30px; } .S0 { margin: 3px 10px 5px 4px; padding: 0px; line-height: 28.8px; min-height: 0px; white-space: pre-wrap; color: rgb(213, 80, 0); font-family: Helvetica, Arial, sans-serif; font-style: normal; font-size: 24px; font-weight: normal; text-align: left;  }
.S1 { margin: 20px 10px 5px 4px; padding: 0px; line-height: 20px; min-height: 0px; white-space: pre-wrap; color: rgb(60, 60, 60); font-family: Helvetica, Arial, sans-serif; font-style: normal; font-size: 20px; font-weight: bold; text-align: left;  }
.S2 { margin: 2px 10px 9px 4px; padding: 0px; line-height: 21px; min-height: 0px; white-space: pre-wrap; color: rgb(0, 0, 0); font-family: Helvetica, Arial, sans-serif; font-style: normal; font-size: 14px; font-weight: normal; text-align: left;  }
.S3 { margin: 2px 10px 9px 4px; padding: 0px; line-height: 21px; min-height: 0px; white-space: pre-wrap; color: rgb(0, 0, 0); font-family: Helvetica, Arial, sans-serif; font-style: normal; font-size: 14px; font-weight: normal; text-align: center;  }
.S4 { margin: 10px 0px 20px; padding-left: 0px; font-family: Helvetica, Arial, sans-serif; font-size: 14px;  }
.S5 { margin-left: 56px; line-height: 21px; min-height: 0px; text-align: left; white-space: pre-wrap;  }
.S6 { margin: 3px 10px 5px 4px; padding: 0px; line-height: 20px; min-height: 0px; white-space: pre-wrap; color: rgb(60, 60, 60); font-family: Helvetica, Arial, sans-serif; font-style: normal; font-size: 20px; font-weight: bold; text-align: left;  }
.CodeBlock { background-color: #F7F7F7; margin: 10px 0 10px 0;}
.S7 { border-left: 1px solid rgb(233, 233, 233); border-right: 1px solid rgb(233, 233, 233); border-top: 1px solid rgb(233, 233, 233); border-bottom: 0px none rgb(0, 0, 0); border-radius: 4px 4px 0px 0px; padding: 6px 45px 0px 13px; line-height: 17.234px; min-height: 18px; white-space: nowrap; color: rgb(0, 0, 0); font-family: Menlo, Monaco, Consolas, "Courier New", monospace; font-size: 14px;  }
.S8 { border-left: 1px solid rgb(233, 233, 233); border-right: 1px solid rgb(233, 233, 233); border-top: 0px none rgb(0, 0, 0); border-bottom: 0px none rgb(0, 0, 0); border-radius: 0px; padding: 0px 45px 0px 13px; line-height: 17.234px; min-height: 18px; white-space: nowrap; color: rgb(0, 0, 0); font-family: Menlo, Monaco, Consolas, "Courier New", monospace; font-size: 14px;  }
.S9 { border-left: 1px solid rgb(233, 233, 233); border-right: 1px solid rgb(233, 233, 233); border-top: 0px none rgb(0, 0, 0); border-bottom: 1px solid rgb(233, 233, 233); border-radius: 0px 0px 4px 4px; padding: 0px 45px 4px 13px; line-height: 17.234px; min-height: 18px; white-space: nowrap; color: rgb(0, 0, 0); font-family: Menlo, Monaco, Consolas, "Courier New", monospace; font-size: 14px;  }
.S10 { margin: 10px 10px 9px 4px; padding: 0px; line-height: 21px; min-height: 0px; white-space: pre-wrap; color: rgb(0, 0, 0); font-family: Helvetica, Arial, sans-serif; font-style: normal; font-size: 14px; font-weight: normal; text-align: left;  }
.S11 { border-left: 1px solid rgb(233, 233, 233); border-right: 1px solid rgb(233, 233, 233); border-top: 1px solid rgb(233, 233, 233); border-bottom: 1px solid rgb(233, 233, 233); border-radius: 4px; padding: 6px 45px 4px 13px; line-height: 17.234px; min-height: 18px; white-space: nowrap; color: rgb(0, 0, 0); font-family: Menlo, Monaco, Consolas, "Courier New", monospace; font-size: 14px;  }
.S12 { border-left: 1px solid rgb(233, 233, 233); border-right: 1px solid rgb(233, 233, 233); border-top: 1px solid rgb(233, 233, 233); border-bottom: 1px solid rgb(233, 233, 233); border-radius: 4px 4px 0px 0px; padding: 6px 45px 4px 13px; line-height: 17.234px; min-height: 18px; white-space: nowrap; color: rgb(0, 0, 0); font-family: Menlo, Monaco, Consolas, "Courier New", monospace; font-size: 14px;  }
.S13 { color: rgb(64, 64, 64); padding: 10px 0px 6px 17px; background: rgb(255, 255, 255) none repeat scroll 0% 0% / auto padding-box border-box; font-family: Menlo, Monaco, Consolas, "Courier New", monospace; font-size: 14px; overflow-x: hidden; line-height: 17.234px;  }
/* Styling that is common to warnings and errors is in diagnosticOutput.css */.embeddedOutputsErrorElement {    min-height: 18px;    max-height: 250px;    overflow: auto;}
.embeddedOutputsErrorElement.inlineElement {}
.embeddedOutputsErrorElement.rightPaneElement {}
/* Styling that is common to warnings and errors is in diagnosticOutput.css */.embeddedOutputsWarningElement{    min-height: 18px;    max-height: 250px;    overflow: auto;}
.embeddedOutputsWarningElement.inlineElement {}
.embeddedOutputsWarningElement.rightPaneElement {}
/* Copyright 2015-2019 The MathWorks, Inc. *//* In this file, styles are not scoped to rtcContainer since they could be in the Dojo Tooltip */.diagnosticMessage-wrapper {    font-family: Menlo, Monaco, Consolas, "Courier New", monospace;    font-size: 12px;}
.diagnosticMessage-wrapper.diagnosticMessage-warningType {    color: rgb(255,100,0);}
.diagnosticMessage-wrapper.diagnosticMessage-warningType a {    color: rgb(255,100,0);    text-decoration: underline;}
.diagnosticMessage-wrapper.diagnosticMessage-errorType {    color: rgb(230,0,0);}
.diagnosticMessage-wrapper.diagnosticMessage-errorType a {    color: rgb(230,0,0);    text-decoration: underline;}
.diagnosticMessage-wrapper .diagnosticMessage-messagePart,.diagnosticMessage-wrapper .diagnosticMessage-causePart {    white-space: pre-wrap;}
.diagnosticMessage-wrapper .diagnosticMessage-stackPart {    white-space: pre;}
.embeddedOutputsTextElement,.embeddedOutputsVariableStringElement {    white-space: pre;    word-wrap:  initial;    min-height: 18px;    max-height: 250px;    overflow: auto;}
.textElement,.rtcDataTipElement .textElement {    padding-top: 3px;}
.embeddedOutputsTextElement.inlineElement,.embeddedOutputsVariableStringElement.inlineElement {}
.inlineElement .textElement {}
.embeddedOutputsTextElement.rightPaneElement,.embeddedOutputsVariableStringElement.rightPaneElement {    min-height: 16px;}
.rightPaneElement .textElement {    padding-top: 2px;    padding-left: 9px;}
.S14 { border-left: 1px solid rgb(233, 233, 233); border-right: 1px solid rgb(233, 233, 233); border-top: 1px solid rgb(233, 233, 233); border-bottom: 0px none rgb(0, 0, 0); border-radius: 0px; padding: 6px 45px 0px 13px; line-height: 17.234px; min-height: 18px; white-space: nowrap; color: rgb(0, 0, 0); font-family: Menlo, Monaco, Consolas, "Courier New", monospace; font-size: 14px;  }
.S15 { border-left: 1px solid rgb(233, 233, 233); border-right: 1px solid rgb(233, 233, 233); border-top: 0px none rgb(0, 0, 0); border-bottom: 1px solid rgb(233, 233, 233); border-radius: 0px; padding: 0px 45px 4px 13px; line-height: 17.234px; min-height: 18px; white-space: nowrap; color: rgb(0, 0, 0); font-family: Menlo, Monaco, Consolas, "Courier New", monospace; font-size: 14px;  }</style></head><body><div class = rtcContent><h1  class = 'S0'><span>Flux Balance Analysis (FBA) and its variants</span></h1><h2  class = 'S1'><span>Author(s): </span><span style=' font-weight: bold;'>Ronan M.T. Fleming,</span><span> Leiden University;</span><span style=' font-weight: bold;'>Vanja Vlasov, LCSB, University of Luxembourg; Thomas Pfau, Systems Biology Group, LSRU, University of Luxembourg</span><span> </span></h2><h2  class = 'S1'><span>Reviewer(s): Ines Thiele,</span><span style=' font-weight: bold;'> Catherine Clancy, National University of Ireland, Galway. Thomas Pfau, Systems Biology Group, LSRU, University of Luxembourg.</span><span> </span></h2><h2  class = 'S1'><span>INTRODUCTION</span></h2><div  class = 'S2'><span>Flux balance analysis (FBA) evaluates the metabolic flux distribution</span><span texencoding="$^1" style="vertical-align:-5px"><img src="" width="8" height="19" /></span><span>, and is one of the most used modelling approaches for metabolic systems. </span></div><div  class = 'S2'><span>The applications of FBA for molecular systems biology include prediction of the growth rates, uptake rates, knockout lethality and product secretion. In FBA, the solution space is constrained by the assumption of a steady-state, under which each internal metabolite is consumed at the same rate as it is produced.</span></div><div  class = 'S2'><span>For the quantitative estimation of the metabolic fluxes, linear programming (LP) can be used to solve the stoichiometric matrix for a given objective function under different constraints. The constraints of the problem depict the space of all eligible possibilities from which an optimal solution can be selected; </span></div><div  class = 'S3'><span texencoding="\begin{array}{ll}
\min\limits _{v} &amp; \ c^{T}v\\
\text{s.t.} &amp; Sv=b,\\
 &amp; l\leq v\leq u,
\end{array}" style="vertical-align:-30px"><img src="" width="97.5" height="72" /></span></div><div  class = 'S2'><span>Equation 1: Formula of standard FBA.</span></div><div  class = 'S2'><span>where </span><span texencoding="$c\in\Re^{n}$" style="vertical-align:-5px"><img src="" width="46.5" height="19" /></span><span> is a parameter vector that linearly combines one or more reaction fluxes to form what is termed the objective function,  and where a </span><span texencoding="$b_{i}&lt;0$" style="vertical-align:-6px"><img src="" width="39.5" height="20" /></span><span>, or  </span><span texencoding="$b_{i}&gt;0$" style="vertical-align:-6px"><img src="" width="39.5" height="20" /></span><span>, represents some fixed output, or input, of the ith molecular species. </span><span texencoding="$S\in\Re^{m\times n}$" style="vertical-align:-5px"><img src="" width="62.5" height="19" /></span><span> is a stoichiometric matrix for </span><span style="font-family: STIXGeneral, STIXGeneral-webfont, serif; font-style: italic; font-weight: normal; color: rgb(0, 0, 0);">m</span><span> molecular species and </span><span style="font-family: STIXGeneral, STIXGeneral-webfont, serif; font-style: italic; font-weight: normal; color: rgb(0, 0, 0);">n</span><span> reactions, and </span><span style="font-family: STIXGeneral, STIXGeneral-webfont, serif; font-style: italic; font-weight: normal; color: rgb(0, 0, 0);">b</span><span> is a vector of known metabolic exchanges. The output of FBA is a particular flux distribution, </span><span style="font-family: STIXGeneral, STIXGeneral-webfont, serif; font-style: italic; font-weight: normal; color: rgb(0, 0, 0);">v</span><span> , which maximises or minimises the objective function and stands between upper and lower bounds, </span><span style="font-family: STIXGeneral, STIXGeneral-webfont, serif; font-style: italic; font-weight: normal; color: rgb(0, 0, 0);">u</span><span> and </span><span style="font-family: STIXGeneral, STIXGeneral-webfont, serif; font-style: italic; font-weight: normal; color: rgb(0, 0, 0);">l</span><span>, respectively.</span></div><div  class = 'S2'><span>There are multiple different variants of FBA, some of which will be demonstrated here:</span></div><ol  class = 'S4'><li  class = 'S5'><span style=' font-weight: bold;'>Standard FBA</span></li><li  class = 'S5'><span style=' font-weight: bold;'>Regularised FBA</span></li><li  class = 'S5'><span style=' font-weight: bold;'>Sparse FBA</span></li><li  class = 'S5'><span style=' font-weight: bold;'>Metabolite dilution FBA (mdFBA)</span></li><li  class = 'S5'><span style=' font-weight: bold;'>Geometric FBA</span></li><li  class = 'S5'><span style=' font-weight: bold;'>Parsimonious enzyme usage Flux Balance Analysis (pFBA)</span></li><li  class = 'S5'><span style=' font-weight: bold;'>Dynamic FBA</span></li><li  class = 'S5'><span style=' font-weight: bold;'>Relaxed FBA</span></li><li  class = 'S5'><span style=' font-weight: bold;'>Flux enrichment analysis (FEA)</span></li></ol><h2  class = 'S6'><span>EQUIPMENT SETUP</span></h2><h2  class = 'S1'><span>Initialise The Cobra Toolbox and set the solver.</span></h2><div  class = 'S2'><span>If necessary, initialise the cobra toolbox:</span></div><div class="CodeBlock"><div class="inlineWrapper"><div  class = 'S7'><span style="white-space: pre"><span > </span><span style="color: rgb(14, 0, 255);">if </span><span >0</span></span></div></div><div class="inlineWrapper"><div  class = 'S8'><span style="white-space: pre"><span >     initCobraToolbox(false) </span><span style="color: rgb(2, 128, 9);">% false, as we don't want to update</span></span></div></div><div class="inlineWrapper"><div  class = 'S9'><span style="white-space: pre"><span > </span><span style="color: rgb(14, 0, 255);">end</span></span></div></div></div><div  class = 'S10'><span>For solving LP problems in a FBA analysis, certain solvers are required and can be set using the </span><span style=' font-family: monospace;'>changeCobraSolver</span><span> function:</span></div><div class="CodeBlock"><div class="inlineWrapper"><div  class = 'S11'><span style="white-space: pre"><span style="color: rgb(2, 128, 9);">% solverOK = changeCobraSolver(solverName, solverType, printLevel, unchecked)</span></span></div></div></div><div  class = 'S10'><span>The present tutoria</span><span>l can ru</span><span>n </span><span>with the GLPK package, which does not require additional installation and configuration. Although, for the analysis of large models is recommended to use the GUROBI package.</span></div><div  class = 'S2'><span>Setup the appropriate solver for the machine you are using by removing the "%" (comment) sign for only the desired solver.</span></div><div class="CodeBlock"><div class="inlineWrapper outputs"><div  class = 'S12'><span style="white-space: pre"><span > changeCobraSolver(</span><span style="color: rgb(170, 4, 249);">'glpk'</span><span >,</span><span style="color: rgb(170, 4, 249);">'all'</span><span >);</span></span></div><div  class = 'S13'><div class="inlineElement eoOutputWrapper embeddedOutputsTextElement scrollableOutput" uid="0857C796" data-testid="output_0" data-width="420" data-height="101" data-hashorizontaloverflow="true" style="width: 450px; max-height: 261px; white-space: pre; font-style: normal; color: rgb(64, 64, 64); font-size: 12px;"><div class="textElement" style="white-space: pre; font-style: normal; color: rgb(64, 64, 64); font-size: 12px;"> &gt; The solver compatibility is not tested with MATLAB R2019b.
 &gt; Solver for LP problems has been set to glpk.
 &gt; The solver compatibility is not tested with MATLAB R2019b.
 &gt; Solver for MILP problems has been set to glpk.
 &gt; Solver glpk not supported for problems of type MIQP. No solver set for this problemtype 
 &gt; Solver glpk not supported for problems of type NLP. No solver set for this problemtype 
 &gt; Solver glpk not supported for problems of type QP. Currently used: qpng </div></div></div></div><div class="inlineWrapper"><div  class = 'S14'><span style="white-space: pre"><span style="color: rgb(2, 128, 9);">% changeCobraSolver('tomlab_cplex','all');</span></span></div></div><div class="inlineWrapper"><div  class = 'S8'><span style="white-space: pre"><span style="color: rgb(2, 128, 9);">% changeCobraASolver('ibm_cplex','all');</span></span></div></div><div class="inlineWrapper"><div  class = 'S9'><span style="white-space: pre"><span style="color: rgb(2, 128, 9);">% changeCobraSolver ('gurobi', 'all');</span></span></div></div></div><h2  class = 'S1'><span>Model Setup</span></h2><div  class = 'S2'><span>Thi</span><span>s tutorial will use the generic model of the human cellular metabolism</span><span texencoding="$ ^2" style="vertical-align:-5px"><img src="" width="8" height="19" /></span><span>, Recon 2.0. Other COBRA models, including </span><span>Recon 3, may also be run with this tutorial. For information on metabolites structures and reactions, and to download the latest COBRA model releases, visit the Virtual Metabolic Human database (VMH, </span><a href = "http://vmh.life)."><span>http://vmh.life).</span></a></div><div  class = 'S2'><span>Before proceeding with the simulations, load the model into the workspace:</span></div><div class="CodeBlock"><div class="inlineWrapper"><div  class = 'S7'><span style="white-space: pre"><span style="color: rgb(14, 0, 255);">global </span><span >CBTDIR</span></span></div></div><div class="inlineWrapper"><div  class = 'S8'><span style="white-space: pre"><span style="color: rgb(2, 128, 9);">%modelFileName = 'Recon2.0model.mat';</span></span></div></div><div class="inlineWrapper"><div  class = 'S8'><span style="white-space: pre"><span >modelFileName = </span><span style="color: rgb(170, 4, 249);">'Recon3DModel_301.mat'</span><span >;</span></span></div></div><div class="inlineWrapper"><div  class = 'S8'><span style="white-space: pre"><span >modelDirectory = getDistributedModelFolder(modelFileName); </span><span style="color: rgb(2, 128, 9);">%Look up the folder for the distributed Models.</span></span></div></div><div class="inlineWrapper"><div  class = 'S8'><span style="white-space: pre"><span >modelFileName= [modelDirectory filesep modelFileName]; </span><span style="color: rgb(2, 128, 9);">% Get the full path. Necessary to be sure, that the right model is loaded</span></span></div></div><div class="inlineWrapper"><div  class = 'S8'><span style="white-space: pre"><span style="color: rgb(14, 0, 255);">if </span><span >0</span></span></div></div><div class="inlineWrapper"><div  class = 'S8'><span style="white-space: pre"><span >    model = readCbModel(modelFileName);</span></span></div></div><div class="inlineWrapper"><div  class = 'S8'><span style="white-space: pre"><span style="color: rgb(14, 0, 255);">else</span></span></div></div><div class="inlineWrapper"><div  class = 'S8'><span style="white-space: pre"><span >    load(</span><span style="color: rgb(170, 4, 249);">'Recon3DModel_301.mat'</span><span >)</span></span></div></div><div class="inlineWrapper"><div  class = 'S8'><span style="white-space: pre"><span >    model = Recon3DModel;</span></span></div></div><div class="inlineWrapper"><div  class = 'S9'><span style="white-space: pre"><span style="color: rgb(14, 0, 255);">end</span></span></div></div></div><div  class = 'S10'><span>In this tutorial we assume, that the cellular objectives include energy production or optimisation of uptake rates and by-product secretion for various physiological functions of the human body.</span></div><h2  class = 'S6'><span>PROCEDURE</span></h2><h2  class = 'S1'><span style=' text-decoration: underline;'>1. Standard FBA</span></h2><div  class = 'S2'><span>Standard FBA predicts an optimal solution for a cellular objective within a given set of constraints on a metabolic network (see Equation 1). Constraints on the network are set by assigning limits on the uptake, consumption or production of metabolites in reactions.</span></div><div  class = 'S2'><span style=' font-weight: bold;'>Timing:</span></div><div  class = 'S2'><span>The time to determine a FBA solution depends on the size of the genome-scale model and is commonly less than a second for a medium sized model.</span></div><div  class = 'S2'><span style=' font-weight: bold;'>Calculating maximal ATP energy production under aerobic conditions:</span></div><div  class = 'S2'><span>For each new simulation, the original model will be copied to a new variable. This preserves the constraints of the original model to perform further simulations with new constraints. Additionally, this method of renaming the model avoids confusion while performing multiple simulations at the same time.</span></div><div class="CodeBlock"><div class="inlineWrapper"><div  class = 'S11'><span style="white-space: pre"><span >modelaerobic = model;</span></span></div></div></div><div  class = 'S10'><span>The ATP demand reacti</span><span>on, i.e., </span><a href = "http://vmh.life/#human/all/DM_atp_c_"><span style=' font-family: monospace;'>DM_atp_c_</span></a><span style=' font-family: monospace;'> </span><span>within the model is a reaction that involves hydrolysis of ATP to ADP, Pi and proton in the cytosol. </span></div><div class="CodeBlock"><div class="inlineWrapper outputs"><div  class = 'S12'><span style="white-space: pre"><span > printRxnFormula(model, </span><span style="color: rgb(170, 4, 249);">'DM_atp_c_'</span><span >);</span></span></div><div  class = 'S13'><div class="inlineElement eoOutputWrapper embeddedOutputsTextElement scrollableOutput" uid="E19D33B1" data-testid="output_1" data-width="420" data-height="18" data-hashorizontaloverflow="true" style="width: 450px; max-height: 261px; white-space: pre; font-style: normal; color: rgb(64, 64, 64); font-size: 12px;"><div class="textElement" style="white-space: pre; font-style: normal; color: rgb(64, 64, 64); font-size: 12px;">DM_atp_c_	h2o[c] + atp[c] 	-&gt;	h[c] + adp[c] + pi[c] </div></div></div></div></div><div  class = 'S10'><span>We will set this reaction</span><span> as our obj</span><span>ec</span><span>tive with the </span><span style=' font-family: monospace;'>'changeObjective'</span><span> command. Maximising the flux through the ATP demand reaction will result in the network producing a maximal amount of ATP (up to the limit of the reaction).</span></div><div class="CodeBlock"><div class="inlineWrapper"><div  class = 'S11'><span style="white-space: pre"><span >modelaerobic = changeObjective (modelaerobic, </span><span style="color: rgb(170, 4, 249);">'DM_atp_c_'</span><span >);</span></span></div></div></div><div  class = 'S10'><span>The glucose and oxygen, in this case, are provid</span><span>ed in high </span><span>am</span><span>ounts for</span><span> ca</span><span>lc</span><span>ulating the flux through ATP demand. </span></div><div  class = 'S2'><span>The </span><span style=' font-family: monospace;'>'changeRxnBounds' </span><span>function changes the flux constraints of the lower ('</span><span style=' font-family: monospace;'>l</span><span>'), upper ('</span><span style=' font-family: monospace;'>u</span><span>'), or both the bounds ('</span><span style=' font-family: monospace;'>b</span><span>'), of the specified reaction. Here, we will change the maximal uptake of glucose to 20 </span><span style="font-family: STIXGeneral, STIXGeneral-webfont, serif; font-style: italic; font-weight: normal; color: rgb(0, 0, 0);">μ</span><span>mol/min/gDW and of oxygen to 1000 </span><span style="font-family: STIXGeneral, STIXGeneral-webfont, serif; font-style: italic; font-weight: normal; color: rgb(0, 0, 0);">μ</span><span>mol/min/gDW. The uptake of oxygen is effectively unconstrainted (i.e. infinity). </span></div><div class="CodeBlock"><div class="inlineWrapper"><div  class = 'S7'><span style="white-space: pre"><span style="color: rgb(2, 128, 9);">% modelaerobic = changeRxnBounds (modelaerobic, 'EX_glc_D[e]', -20, 'l'); % For Recon 3.0 uncomment these lines and </span></span></div></div><div class="inlineWrapper"><div  class = 'S8'><span style="white-space: pre"><span style="color: rgb(2, 128, 9);">%modelaerobic = changeRxnBounds (modelaerobic, 'EX_o2[e]', -1000, 'l'); % comment the lines below.</span></span></div></div><div class="inlineWrapper"><div  class = 'S8'><span style="white-space: pre"><span >modelaerobic = changeRxnBounds (modelaerobic, </span><span style="color: rgb(170, 4, 249);">'EX_glc_D[e]'</span><span >, -20, </span><span style="color: rgb(170, 4, 249);">'l'</span><span >);  </span></span></div></div><div class="inlineWrapper"><div  class = 'S9'><span style="white-space: pre"><span >modelaerobic = changeRxnBounds (modelaerobic, </span><span style="color: rgb(170, 4, 249);">'EX_o2[e]'</span><span >, -1000, </span><span style="color: rgb(170, 4, 249);">'l'</span><span >); </span></span></div></div></div><div  class = 'S10'><span>The function </span><span style=' font-family: monospace;'>optimizeCbModel</span><span> calculates one of the op</span><span>timal solu</span><span>tions for a (maximum or minimum) objective reaction wit</span><span>hin</span><span> t</span><span>he defined solution space. In the above example, the maximal flux through the </span><a href = "http://vmh.life/#human/all/DM_atp_c_"><span style=' font-family: monospace;'>DM_atp_c_</span></a><span> is desired. </span></div><div class="CodeBlock"><div class="inlineWrapper outputs"><div  class = 'S12'><span style="white-space: pre"><span >FBAaerobic = optimizeCbModel (modelaerobic, </span><span style="color: rgb(170, 4, 249);">'max'</span><span >)</span></span></div><div  class = 'S13'><div class="inlineElement eoOutputWrapper embeddedOutputsVariableStringElement" uid="976F7EB0" data-testid="output_2" data-width="420" data-height="258" data-hashorizontaloverflow="false" style="width: 450px; max-height: 269px; white-space: pre; font-style: normal; color: rgb(64, 64, 64); font-size: 12px;"><div class="textElement" style="white-space: pre; font-style: normal; color: rgb(64, 64, 64); font-size: 12px;"><div style="white-space: pre; font-style: normal; color: rgb(64, 64, 64); font-size: 12px;"><span class="variableNameElement" style="white-space: pre; font-style: normal; color: rgb(64, 64, 64); font-size: 12px;">FBAaerobic = <span class="headerElement" style="white-space: pre; font-style: italic; color: rgb(179, 179, 179); font-size: 12px;">struct with fields:</span></span></div><div style="white-space: pre; font-style: normal; color: rgb(64, 64, 64); font-size: 12px;">         full: [10600×1 double]
          obj: 1000
        rcost: [10600×1 double]
         dual: [5835×1 double]
        slack: [5835×1 double]
       solver: 'glpk'
    algorithm: 'default'
         stat: 1
     origStat: 5
         time: 0.7090
        basis: []
            f: 1000
            x: [10600×1 double]
            v: [10600×1 double]
            w: [10600×1 double]
            y: [5835×1 double]
            s: [5835×1 double]
</div></div></div></div></div></div><ul  class = 'S4'><li  class = 'S5'><span style=' font-weight: bold;'>Anticipated results</span></li></ul><div  class = 'S2'><span>When oxygen and all carbon sources (internal and external) are provided the flux through ATP demand reaction can reach its maximum rate of 1000 </span><span style="font-family: STIXGeneral, STIXGeneral-webfont, serif; font-style: italic; font-weight: normal; color: rgb(0, 0, 0);">μ</span><span>mol/min/gDW. </span></div><ul  class = 'S4'><li  class = 'S5'><span style=' font-weight: bold;'>Troubleshooting</span></li></ul><div  class = 'S2'><span>If there are multiple carbon sources available in the model, it may be necessary to specify more constraints in order to examine the effect of a single carbon source on ATP production.  </span></div><div  class = 'S2'><span>To avoid this issue, all external carbon sources need to be closed with the exception of the single carbon source of interest.</span></div><div class="CodeBlock"><div class="inlineWrapper"><div  class = 'S7'><span style="white-space: pre"><span style="color: rgb(2, 128, 9);">%Closing the uptake of all energy and oxygen sources</span></span></div></div><div class="inlineWrapper"><div  class = 'S8'><span style="white-space: pre"><span >[exchBool,uptBool] = findExcRxns(model);</span></span></div></div><div class="inlineWrapper"><div  class = 'S8'><span style="white-space: pre"><span >uptakes = model.rxns(uptBool);</span></span></div></div><div class="inlineWrapper"><div  class = 'S8'></div></div><div class="inlineWrapper"><div  class = 'S8'><span style="white-space: pre"><span style="color: rgb(2, 128, 9);">% If you use Recon3.0 model, then:</span></span></div></div><div class="inlineWrapper"><div  class = 'S8'><span style="white-space: pre"><span style="color: rgb(2, 128, 9);">% modelalter = model;</span></span></div></div><div class="inlineWrapper"><div  class = 'S8'><span style="white-space: pre"><span style="color: rgb(2, 128, 9);">% modelalter = changeRxnBounds(modelalter, uptakes, 0, 'b');</span></span></div></div><div class="inlineWrapper"><div  class = 'S8'><span style="white-space: pre"><span style="color: rgb(2, 128, 9);">% modelalter = changeRxnBounds(modelalter, 'EX_HC00250[e]', -1000, 'l');</span></span></div></div><div class="inlineWrapper"><div  class = 'S8'></div></div><div class="inlineWrapper"><div  class = 'S8'><span style="white-space: pre"><span style="color: rgb(2, 128, 9);">% The alternative way to do that, in case you were using another large model, </span></span></div></div><div class="inlineWrapper"><div  class = 'S8'><span style="white-space: pre"><span style="color: rgb(2, 128, 9);">% that does not contain defined Subsystem is</span></span></div></div><div class="inlineWrapper"><div  class = 'S8'><span style="white-space: pre"><span style="color: rgb(2, 128, 9);">% to find uptake exchange reactions with following codes:</span></span></div></div><div class="inlineWrapper"><div  class = 'S8'><span style="white-space: pre"><span style="color: rgb(2, 128, 9);">% [selExc, selUpt] = findExcRxns(model);</span></span></div></div><div class="inlineWrapper"><div  class = 'S8'><span style="white-space: pre"><span style="color: rgb(2, 128, 9);">% uptakes1 = model.rxns(selUpt);</span></span></div></div><div class="inlineWrapper"><div  class = 'S8'></div></div><div class="inlineWrapper"><div  class = 'S8'><span style="white-space: pre"><span style="color: rgb(2, 128, 9);">% Selecting from the exchange uptake reactions those </span></span></div></div><div class="inlineWrapper"><div  class = 'S8'><span style="white-space: pre"><span style="color: rgb(2, 128, 9);">% which contain at least 1 carbon in the metabolites included in the reaction:</span></span></div></div><div class="inlineWrapper"><div  class = 'S8'><span style="white-space: pre"><span > subuptakeModel = extractSubNetwork(model, uptakes);</span></span></div></div><div class="inlineWrapper"><div  class = 'S8'><span style="white-space: pre"><span > hiCarbonRxns = findCarbonRxns(subuptakeModel,1);</span></span></div></div><div class="inlineWrapper"><div  class = 'S8'><span style="white-space: pre"><span style="color: rgb(2, 128, 9);">% Closing the uptake of all the carbon sources</span></span></div></div><div class="inlineWrapper"><div  class = 'S8'><span style="white-space: pre"><span > modelalter = model;</span></span></div></div><div class="inlineWrapper"><div  class = 'S8'><span style="white-space: pre"><span > modelalter = changeRxnBounds(modelalter, hiCarbonRxns, 0, </span><span style="color: rgb(170, 4, 249);">'b'</span><span >);</span></span></div></div><div class="inlineWrapper"><div  class = 'S8'><span style="white-space: pre"><span style="color: rgb(2, 128, 9);">% Closing other oxygen and energy sources. Use the following lines for recon2, or uncomment the lines below for recon3</span></span></div></div><div class="inlineWrapper"><div  class = 'S8'><span style="white-space: pre"><span > exoxygen = {</span><span style="color: rgb(170, 4, 249);">'EX_adp[e]'</span><span >,   </span><span style="color: rgb(170, 4, 249);">'EX_amp[e]'</span><span >,    </span><span style="color: rgb(170, 4, 249);">'EX_atp[e]'</span><span >,    </span><span style="color: rgb(170, 4, 249);">'EX_co2[e]'</span><span >,    </span><span style="color: rgb(170, 4, 249);">'EX_coa[e]'</span><span >,    </span><span style="color: rgb(170, 4, 249);">'EX_fad[e]'</span><span >,    </span><span style="color: rgb(170, 4, 249);">'EX_fe2[e]'</span><span >,</span><span style="color: rgb(14, 0, 255);">...</span></span></div></div><div class="inlineWrapper"><div  class = 'S8'><span style="white-space: pre"><span >    </span><span style="color: rgb(170, 4, 249);">'EX_fe3[e]'</span><span >,    </span><span style="color: rgb(170, 4, 249);">'EX_gdp[e]'</span><span >,   </span><span style="color: rgb(170, 4, 249);">'EX_gmp[e]'</span><span >,    </span><span style="color: rgb(170, 4, 249);">'EX_gtp[e]'</span><span >,    </span><span style="color: rgb(170, 4, 249);">'EX_h[e]'</span><span >,    </span><span style="color: rgb(170, 4, 249);">'EX_h2o[e]'</span><span >,    </span><span style="color: rgb(170, 4, 249);">'EX_h2o2[e]'</span><span >,</span><span style="color: rgb(14, 0, 255);">...</span></span></div></div><div class="inlineWrapper"><div  class = 'S8'><span style="white-space: pre"><span >    </span><span style="color: rgb(170, 4, 249);">'EX_nad[e]'</span><span >,    </span><span style="color: rgb(170, 4, 249);">'EX_nadp[e]'</span><span >,    </span><span style="color: rgb(170, 4, 249);">'EX_no[e]'</span><span >,    </span><span style="color: rgb(170, 4, 249);">'EX_no2[e]'</span><span >,    </span><span style="color: rgb(170, 4, 249);">'EX_o2s[e]'</span><span >};</span></span></div></div><div class="inlineWrapper"><div  class = 'S9'><span style="white-space: pre"><span >modelalter = changeRxnBounds (modelalter, exoxygen, 0, </span><span style="color: rgb(170, 4, 249);">'l'</span><span >);</span></span></div></div></div><div  class = 'S10'><span style=' font-weight: bold;'>Calculating maximum ATP energy production under anaerobic and glucose only conditions:</span></div><div class="CodeBlock"><div class="inlineWrapper"><div  class = 'S7'><span style="white-space: pre"><span >modelanaerobic = modelalter;</span></span></div></div><div class="inlineWrapper"><div  class = 'S8'><span style="white-space: pre"><span style="color: rgb(2, 128, 9);">% modelaerobic = changeRxnBounds (modelaerobic, 'EX_glc_D[e]', -20, 'l'); % For Recon 3.0 uncomment these lines and </span></span></div></div><div class="inlineWrapper"><div  class = 'S8'><span style="white-space: pre"><span style="color: rgb(2, 128, 9);">%modelaerobic = changeRxnBounds (modelaerobic, 'EX_o2[e]', -1000, 'l'); % comment the lines below.</span></span></div></div><div class="inlineWrapper"><div  class = 'S8'><span style="white-space: pre"><span >modelanaerobic = changeRxnBounds(modelanaerobic, </span><span style="color: rgb(170, 4, 249);">'EX_glc_D[e]'</span><span >,-20,</span><span style="color: rgb(170, 4, 249);">'l'</span><span >);</span></span></div></div><div class="inlineWrapper"><div  class = 'S8'><span style="white-space: pre"><span >modelanaerobic = changeRxnBounds (modelanaerobic, </span><span style="color: rgb(170, 4, 249);">'EX_o2[e]'</span><span >, 0, </span><span style="color: rgb(170, 4, 249);">'l'</span><span >);</span></span></div></div><div class="inlineWrapper"><div  class = 'S8'><span style="white-space: pre"><span >modelanaerobic = changeObjective(modelanaerobic,</span><span style="color: rgb(170, 4, 249);">'DM_atp_c_'</span><span >);</span></span></div></div><div class="inlineWrapper outputs"><div  class = 'S15'><span style="white-space: pre"><span >FBAanaerob = optimizeCbModel(modelanaerobic,</span><span style="color: rgb(170, 4, 249);">'max'</span><span >)</span></span></div><div  class = 'S13'><div class="inlineElement eoOutputWrapper embeddedOutputsVariableStringElement" uid="5E730D0A" data-testid="output_3" data-width="420" data-height="258" data-hashorizontaloverflow="false" style="width: 450px; max-height: 269px; white-space: pre; font-style: normal; color: rgb(64, 64, 64); font-size: 12px;"><div class="textElement" style="white-space: pre; font-style: normal; color: rgb(64, 64, 64); font-size: 12px;"><div style="white-space: pre; font-style: normal; color: rgb(64, 64, 64); font-size: 12px;"><span class="variableNameElement" style="white-space: pre; font-style: normal; color: rgb(64, 64, 64); font-size: 12px;">FBAanaerob = <span class="headerElement" style="white-space: pre; font-style: italic; color: rgb(179, 179, 179); font-size: 12px;">struct with fields:</span></span></div><div style="white-space: pre; font-style: normal; color: rgb(64, 64, 64); font-size: 12px;">         full: [10600×1 double]
          obj: 1.0000e+03
        rcost: [10600×1 double]
         dual: [5835×1 double]
        slack: [5835×1 double]
       solver: 'glpk'
    algorithm: 'default'
         stat: 1
     origStat: 5
         time: 0.6735
        basis: []
            f: 1.0000e+03
            x: [10600×1 double]
            v: [10600×1 double]
            w: [10600×1 double]
            y: [5835×1 double]
            s: [5835×1 double]
</div></div></div></div></div></div><ul  class = 'S4'><li  class = 'S5'><span style=' font-weight: bold;'>Anticipated results</span></li></ul><div  class = 'S2'><span>Compared to the aerobic condition, anaerobic condition with only glucose as an energy source has reduced flux through ATP demand (82 </span><span style="font-family: STIXGeneral, STIXGeneral-webfont, serif; font-style: italic; font-weight: normal; color: rgb(0, 0, 0);">μ</span><span>mol/min/gDW), signifying the need to oxygen to run the oxidative phosphorylation. The results are dependant on the model you are using. For Recon 3.0, under anaerobic conditions with only glucose as an energy source, the flux for ATP demand is 40 </span><span style="font-family: STIXGeneral, STIXGeneral-webfont, serif; font-style: italic; font-weight: normal; color: rgb(0, 0, 0);">μ</span><span>mol/min/gDW.</span></div><h2  class = 'S6'><span style=' text-decoration: underline;'>2. Quadratically regularised FBA</span></h2><div  class = 'S2'><span>Regularised FBA calculates the optimal solution of a linear objective function, and finds the unique optimal flux vector that minimises the Euclidean norm of the flux vector, that is</span></div><div  class = 'S3'><span texencoding="\begin{array}{ll}
\min\limits _{v} &amp; \ c^{T}v + \eta ||v|| \\
\text{s.t.} &amp; Sv=b,\\
 &amp; l\leq v\leq u,\\
\end{array}" style="vertical-align:-30px"><img src="" width="116" height="72" /></span></div><div  class = 'S2'><span>If the scalar parameter </span><span style="font-family: STIXGeneral, STIXGeneral-webfont, serif; font-style: italic; font-weight: normal; color: rgb(0, 0, 0);">η</span><span> is small, then </span><span texencoding="c^Tv" style="vertical-align:-5px"><img src="" width="23" height="19" /></span><span> dominates and the same optimal objective as in FBA is attained, with the advantage that a unique flux vector is returned - one that minimises the the sum of the squares of the flux vector.</span></div><div  class = 'S2'><span>In practice, for genome-scale, mono-scale models (e.g. not ME models) the following value for the scalar parameter works fine</span></div><div class="CodeBlock"><div class="inlineWrapper"><div  class = 'S11'><span style="white-space: pre"><span >eta = 1e-6;</span></span></div></div></div><div  class = 'S10'><span>This is a quadratic optimisation problem, so a 'QP' solver is required. The COBRA Toolbox v3 comes shipped with the solver PDCO: Primal-Dual interior method for Convex Objectives (further info here: </span><a href = "https://web.stanford.edu/group/SOL/software/pdco/"><span>link</span></a><span>).</span></div><div class="CodeBlock"><div class="inlineWrapper"><div  class = 'S11'><span style="white-space: pre"><span >[solverOK, solverInstalled] = changeCobraSolver(</span><span style="color: rgb(170, 4, 249);">'pdco'</span><span >,</span><span style="color: rgb(170, 4, 249);">'QP'</span><span >);</span></span></div></div></div><div  class = 'S10'><span>Solve a quadratically regularised FBA problem:</span></div><div class="CodeBlock"><div class="inlineWrapper"><div  class = 'S11'><span style="white-space: pre"><span >FBAsolution = optimizeCbModel(model,</span><span style="color: rgb(170, 4, 249);">'max'</span><span >,eta);</span></span></div></div></div><h2  class = 'S1'><span style=' text-decoration: underline;'>2. Sparse FBA</span></h2><div  class = 'S2'><span>Sparse FBA calculates the optimal solution of a linear objective function and finds the smallest set of reactions that can carry flux to achieve the objective. Sparse FBA minimises the number of active reactions by keeping same maximal objective;</span></div><div  class = 'S3'><span texencoding="\begin{array}{ll}
\min\limits _{v} &amp; \ ||v||_0 \\
\text{s.t.} &amp; Sv=b,\\
 &amp; l\leq v\leq u,\\
&amp; c^{T}v=\rho^*
\end{array}" style="vertical-align:-41px"><img src="" width="97.5" height="93" /></span></div><div  class = 'S2'><span>Equation 2: Formula of Sparse FBA.</span></div><div  class = 'S2'><span>where the last constraint is optional and represents the requirement to satisfy an optimal objective value </span><span texencoding="\rho^*" style="vertical-align:-5px"><img src="" width="16.5" height="19" /></span><span> derived from any solution to a FBA problem. This approach is used to check for minimal sets of reactions that either should be active or should not be active in a flux balance model that is representative of a biochemical network.</span></div><div class="CodeBlock"><div class="inlineWrapper"><div  class = 'S7'><span style="white-space: pre"><span style="color: rgb(2, 128, 9);">% [vSparse, sparseRxnBool, essentialRxnBool]  = sparseFBA(model, osenseStr,...</span></span></div></div><div class="inlineWrapper"><div  class = 'S9'><span style="white-space: pre"><span style="color: rgb(2, 128, 9);">%  checkMinimalSet, checkEssentialSet, zeroNormApprox)</span></span></div></div></div><div  class = 'S10'><span>As an optional input, there are different appoximation types of zero-norm (only available when </span><span style=' font-family: monospace;'>minNorm = 'zero'</span><span>). Default is </span><span style=' font-family: monospace;'>cappedL1</span><span>.</span></div><div class="CodeBlock"><div class="inlineWrapper"><div  class = 'S7'><span style="white-space: pre"><span style="color: rgb(2, 128, 9);">% Other types of zero-norm:</span></span></div></div><div class="inlineWrapper"><div  class = 'S8'><span style="white-space: pre"><span style="color: rgb(2, 128, 9);">%  * 'cappedL1' : Capped-L1 norm</span></span></div></div><div class="inlineWrapper"><div  class = 'S8'><span style="white-space: pre"><span style="color: rgb(2, 128, 9);">%  * 'exp'      : Exponential function</span></span></div></div><div class="inlineWrapper"><div  class = 'S8'><span style="white-space: pre"><span style="color: rgb(2, 128, 9);">%  * 'log'      : Logarithmic function</span></span></div></div><div class="inlineWrapper"><div  class = 'S8'><span style="white-space: pre"><span style="color: rgb(2, 128, 9);">%  * 'SCAD'     : SCAD function</span></span></div></div><div class="inlineWrapper"><div  class = 'S8'><span style="white-space: pre"><span style="color: rgb(2, 128, 9);">%  * 'lp-'      : :math:`L_p` norm with :math:`p &lt; 0`</span></span></div></div><div class="inlineWrapper"><div  class = 'S8'><span style="white-space: pre"><span style="color: rgb(2, 128, 9);">%  * 'lp+'      : :math:`L_p` norm with :math:`0 &lt; p &lt; 1`</span></span></div></div><div class="inlineWrapper"><div  class = 'S8'><span style="white-space: pre"><span style="color: rgb(2, 128, 9);">%  * 'l1'       : L1 norm</span></span></div></div><div class="inlineWrapper"><div  class = 'S9'><span style="white-space: pre"><span style="color: rgb(2, 128, 9);">%  * 'all'      : try all approximations and return the best result</span></span></div></div></div><div  class = 'S10'><span style=' font-weight: bold;'>Timing:</span></div><div  class = 'S2'><span>The time to determine a </span><span style=' font-family: monospace;'>sparseFBA()</span><span> solution depends on the size of the genome-scale model and is taking from </span><span mathmlencoding="&lt;math xmlns=&quot;http://www.w3.org/1998/Math/MathML&quot; display=&quot;inline&quot;&gt;&lt;mrow&gt;&lt;mo stretchy=&quot;false&quot;&gt;&amp;lt;&lt;/mo&gt;&lt;mn&gt;1&lt;/mn&gt;&lt;/mrow&gt;&lt;/math&gt;" style="vertical-align:-5px"><img src="" width="25" height="18" /></span><span> second for a 1,000 reaction model, to </span><span mathmlencoding="&lt;math xmlns=&quot;http://www.w3.org/1998/Math/MathML&quot; display=&quot;inline&quot;&gt;&lt;mrow&gt;&lt;mo&gt;&amp;lt;&lt;/mo&gt;&lt;mn&gt;2&lt;/mn&gt;&lt;/mrow&gt;&lt;/math&gt;" style="vertical-align:-5px"><img src="" width="25" height="18" /></span><span> seconds for a model with more than 10,000 reactions.</span></div><div  class = 'S2'><span style=' font-weight: bold;'>Calculating maximal ATP energy production under anaerobic and glucose only conditions:</span></div><div class="CodeBlock"><div class="inlineWrapper"><div  class = 'S7'><span style="white-space: pre"><span >modelspar = modelalter;</span></span></div></div><div class="inlineWrapper"><div  class = 'S8'><span style="white-space: pre"><span style="color: rgb(2, 128, 9);">% For Recon3.0 model</span></span></div></div><div class="inlineWrapper"><div  class = 'S8'><span style="white-space: pre"><span style="color: rgb(2, 128, 9);">% modelspar = changeRxnBounds (modelspar, 'EX_glc_D[e]', -20, 'l');</span></span></div></div><div class="inlineWrapper"><div  class = 'S8'><span style="white-space: pre"><span style="color: rgb(2, 128, 9);">% modelspar = changeRxnBounds (modelspar, 'EX_o2[e]', 0, 'l');</span></span></div></div><div class="inlineWrapper"><div  class = 'S8'><span style="white-space: pre"><span >modelspar = changeRxnBounds(modelspar, </span><span style="color: rgb(170, 4, 249);">'EX_glc_D[e]'</span><span >,-20,</span><span style="color: rgb(170, 4, 249);">'l'</span><span >);</span></span></div></div><div class="inlineWrapper"><div  class = 'S8'><span style="white-space: pre"><span >modelspar = changeRxnBounds (modelspar, </span><span style="color: rgb(170, 4, 249);">'EX_o2[e]'</span><span >, 0, </span><span style="color: rgb(170, 4, 249);">'l'</span><span >);</span></span></div></div><div class="inlineWrapper"><div  class = 'S8'><span style="white-space: pre"><span >modelspar = changeObjective(modelspar, </span><span style="color: rgb(170, 4, 249);">'DM_atp_c_'</span><span >);</span></span></div></div><div class="inlineWrapper outputs"><div  class = 'S15'><span style="white-space: pre"><span >[vSparse, sparseRxnBool, essentialRxnBool] = sparseFBA(modelspar, </span><span style="color: rgb(170, 4, 249);">'max'</span><span >);</span></span></div><div  class = 'S13'><div class="inlineElement eoOutputWrapper embeddedOutputsTextElement" uid="56726BA3" data-testid="output_4" data-width="420" data-height="18" data-hashorizontaloverflow="false" style="width: 450px; max-height: 261px; white-space: pre; font-style: normal; color: rgb(64, 64, 64); font-size: 12px;"><div class="textElement" style="white-space: pre; font-style: normal; color: rgb(64, 64, 64); font-size: 12px;">---FBA---</div></div></div></div></div><ul  class = 'S4'><li  class = 'S5'><span style=' font-weight: bold;'>Anticipated results:</span></li></ul><div  class = 'S2'><span>Commonly, a sparse FBA solution will have much smaller number of active reactions compared to a standard FBA on the same model with same objective function. The outputs </span><span style=' font-family: monospace;'>sparseRxnBool </span><span>and </span><span style=' font-family: monospace;'>essentialRxnBool</span><span> return vectors with 1 and 0's, with sparse and essential reactions respectively.</span></div><div  class = 'S2'><span>Display the sparse flux solution, but only the non-zero fluxes. This is convenient with the </span><span style=' font-family: monospace;'>printFluxVector</span><span> function.</span></div><div class="CodeBlock"><div class="inlineWrapper"><div  class = 'S7'><span style="white-space: pre"><span >nonZeroFlag = 1;</span></span></div></div><div class="inlineWrapper outputs"><div  class = 'S15'><span style="white-space: pre"><span >printFluxVector(model, vSparse, nonZeroFlag)</span></span></div><div  class = 'S13'><div class="inlineElement eoOutputWrapper embeddedOutputsTextElement scrollableOutput" uid="4EAC165E" data-testid="output_5" data-width="420" data-height="1039" data-hashorizontaloverflow="false" style="width: 450px; max-height: 261px; white-space: pre; font-style: normal; color: rgb(64, 64, 64); font-size: 12px;"><div class="textElement" style="white-space: pre; font-style: normal; color: rgb(64, 64, 64); font-size: 12px;">10FTHFtm            	         250
4ABUTtm             	        -250
ABTArm              	        -250
ACONTm              	         250
AKGDm               	         250
ATPtm               	         250
CSm                 	         250
DCK1m               	        -250
FOLt2               	         500
GALt1r              	         250
GLUDxm              	         250
GLYtm               	        -500
GTHDH               	         500
H2CO3Dm             	         250
H2Otm               	         250
HPYRRy              	         500
ICDHxm              	         250
MTHFCm              	        -250
MTHFD2m             	        -250
NH4t3r              	        -750
PCm                 	         250
PDHm                	         250
SUCOASm             	        -250
THFtm               	        -500
TRIOK               	         500
EX_HC00250[e]       	        -500
EX_HC02160[e]       	         250
EX_HC02161[e]       	        -250
r0021               	        -500
r0027               	         250
r0081               	         500
r0160               	         500
r0178               	        -250
r0193               	        -500
r0377               	         250
r0838               	         250
r0940               	        -500
r1411               	         250
r1434               	         500
RE2605C             	         500
RE2898C             	         500
EX_oh1[e]           	        -500
FOLt2le             	        -500
DM_Lcystin          	         250
r0295               	         500
MLTHFtm             	         250
sink_thf[c]         	        1000
DM_atp_c_           	        1000
ENO                 	        1000
EX_h[e]             	         750
EX_h2o[e]           	         250
FBA                 	         250
GAPD                	        1000
GHMT2r              	        -500
MTHFC               	         250
MTHFD               	         250
PGI                 	         250
PGK                 	       -1000
PGM                 	       -1000
PGMT                	         250
TPI                 	         250
UDPG4E              	        -250
EX_nh4[e]           	        -750
GALK                	         250
PFK                 	         250
PYK                 	        1000
r0392               	        -500
UGLT                	         250
MLTHFte3            	        1000
EX_mlthf[e]         	       -1000
HMR_6607            	        -500
HMR_8475            	        -250
DM_4abut[c]         	         250
DCMPtm              	        -250</div></div></div></div></div><h2  class = 'S6'><span style=' text-decoration: underline;'>3. Metabolite dilution flux balance analysis (mdFBA)</span></h2><div  class = 'S2'><span>This is a variant of FBA for predicting metabolic flux distributions by accounting for growth-associated dilution of all metabolites in a context-dependent manner</span><span texencoding="$^3" style="vertical-align:-5px"><img src="" width="8" height="19" /></span><span>.</span></div><div  class = 'S2'><span>A solution from the function </span><span style=' font-family: monospace;'>mdFBA</span><span> supports that all metabolites used in any reaction of the solution can either be produced by the network or taken up from the surrounding medium.</span></div><div  class = 'S2'><span style=' font-weight: bold;'>Timing:</span></div><div  class = 'S2'><span>Since this is a MIXED Integer Problem it can take a long time to solve.</span></div><div  class = 'S2'><span style=' font-weight: bold;'>Calculating ATP energy production under aerobic condition using mdFBA:</span></div><div  class = 'S2'><span>In this function, there is an optional output </span><span style=' font-family: monospace;'>newActives</span><span>, th</span><span>at represent reactions that are only active in this analysis.</span></div><div class="CodeBlock"><div class="inlineWrapper"><div  class = 'S7'><span style="white-space: pre"><span style="color: rgb(2, 128, 9);">% The valid solution can be produced with the Recon 3.0 model</span></span></div></div><div class="inlineWrapper"><div  class = 'S8'><span style="white-space: pre"><span style="color: rgb(2, 128, 9);">% modelmd = model;</span></span></div></div><div class="inlineWrapper"><div  class = 'S8'><span style="white-space: pre"><span style="color: rgb(2, 128, 9);">% modelmd = changeRxnBounds(modelmd, 'EX_glc_D[e]',-20,'l');</span></span></div></div><div class="inlineWrapper"><div  class = 'S8'><span style="white-space: pre"><span style="color: rgb(2, 128, 9);">% modelmd = changeRxnBounds (modelmd, 'EX_o2[e]', -1000, 'l');</span></span></div></div><div class="inlineWrapper"><div  class = 'S8'><span style="white-space: pre"><span style="color: rgb(2, 128, 9);">% modelmd = changeObjective(modelmd, 'DM_atp_c_');</span></span></div></div><div class="inlineWrapper"><div  class = 'S8'></div></div><div class="inlineWrapper"><div  class = 'S9'><span style="white-space: pre"><span style="color: rgb(2, 128, 9);">% [sol, newActives] = mdFBA(modelmd)</span></span></div></div></div><ul  class = 'S4'><li  class = 'S5'><span style=' font-weight: bold;'>Troubleshooting:</span></li></ul><div  class = 'S2'><span>When a model do</span><span>es not have a feasible solution, add the input:  </span><span style=' font-family: monospace;'>'getInvalidSolution', true</span><span>.</span></div><div class="CodeBlock"><div class="inlineWrapper"><div  class = 'S7'><span style="white-space: pre"><span style="color: rgb(2, 128, 9);">% clear modelmd</span></span></div></div><div class="inlineWrapper"><div  class = 'S8'><span style="white-space: pre"><span >modelnosol = modelalter;</span></span></div></div><div class="inlineWrapper"><div  class = 'S8'><span style="white-space: pre"><span >modelnosol = changeObjective(modelnosol, </span><span style="color: rgb(170, 4, 249);">'DM_atp_c_'</span><span >);</span></span></div></div><div class="inlineWrapper outputs"><div  class = 'S15'><span style="white-space: pre"><span >[sol, newActives] = mdFBA(modelnosol,  </span><span style="color: rgb(170, 4, 249);">'getInvalidSolution'</span><span >, true)</span></span></div><div  class = 'S13'><div class="inlineElement eoOutputWrapper embeddedOutputsTextElement" uid="26869E5A" data-testid="output_6" data-width="420" data-height="115" data-hashorizontaloverflow="false" style="width: 450px; max-height: 261px; white-space: pre; font-style: normal; color: rgb(64, 64, 64); font-size: 12px;"><div class="textElement" style="white-space: pre; font-style: normal; color: rgb(64, 64, 64); font-size: 12px;">sol =

     []


newActives =

  0×0 empty cell array</div></div></div></div></div><div  class = 'S2'><span>Sometimes when an FBA of a model with the same objective function and constraints is run many times, or using different LP logarithm, we may get different set of solutions for individual reactions. In other words, there are different sets of </span><span style=' font-family: monospace;'>'FBAsolution.v'</span><span> values (fluxes of the reactions) and still get the same objective function value </span><span style=' font-family: monospace;'>'f'</span><span>. That is, the opitmal solution vector is not unique. This can create difficulty when investigating the changes to fluxes between two different conditions. In this case a unique solution is required to compare the changes to fluxes.</span></div><div  class = 'S2'><span>This issue can be solved with either of the following the methods</span><span style=' font-family: monospace;'> </span></div><ul  class = 'S4'><li  class = 'S5'><span style=' font-family: monospace;'>geometricFBA</span><span>, which provides a standard, central and reproducible solution, or</span></li><li  class = 'S5'><span style=' font-family: monospace;'>pFBA</span><span>, which provides a solution based on the minimal fluxes through the model, and classify each gene according to how it contributes to the optimal solution.</span></li></ul><h2  class = 'S1'><span style=' text-decoration: underline;'>4. Geometric FBA</span></h2><div  class = 'S2'><span>The geometric FBA solves the smallest frame that contains all sets of optimal FBA solutions and posts a set of multiple linear programming problems</span><span texencoding="$^4" style="vertical-align:-5px"><img src="" width="8" height="19" /></span><span>.</span></div><div  class = 'S2'><span>This FBA analysis applies iterations, where by each iteration reduces the permissible solution space. After a finite number of iterations, it resolves one single solution of the flux distribution.</span></div><div class="CodeBlock"><div class="inlineWrapper"><div  class = 'S7'><span style="white-space: pre"><span style="color: rgb(2, 128, 9);">% USAGE:</span></span></div></div><div class="inlineWrapper"><div  class = 'S9'><span style="white-space: pre"><span style="color: rgb(2, 128, 9);">% flux = geometricFBA(model, varargin)</span></span></div></div></div><div  class = 'S10'><span style=' font-weight: bold;'>Timing:</span></div><div  class = 'S2'><span>The time to determine </span><span>a geometric FBA solution depends on the size of the genome-scale model and the number of iterations. For a model with more than 10,000 reactions and several iterations takes </span><span mathmlencoding="&lt;math xmlns=&quot;http://www.w3.org/1998/Math/MathML&quot; display=&quot;inline&quot;&gt;&lt;mrow&gt;&lt;mo&gt;&amp;geq;&lt;/mo&gt;&lt;mn&gt;30&lt;/mn&gt;&lt;/mrow&gt;&lt;/math&gt;" style="vertical-align:-5px"><img src="" width="32.5" height="18" /></span><span>minutes.</span></div><div  class = 'S2'><span style=' font-weight: bold;'>Calculating ATP energy production under anaerobic conditions using geometric FBA:</span></div><div class="CodeBlock"><div class="inlineWrapper"><div  class = 'S7'><span style="white-space: pre"><span >modelgeo = modelalter;</span></span></div></div><div class="inlineWrapper"><div  class = 'S8'><span style="white-space: pre"><span style="color: rgb(2, 128, 9);">% For Recon3.0 model</span></span></div></div><div class="inlineWrapper"><div  class = 'S8'><span style="white-space: pre"><span style="color: rgb(2, 128, 9);">% modelgeo = changeRxnBounds (modelgeo, 'EX_glc_D[e]', -20, 'l');</span></span></div></div><div class="inlineWrapper"><div  class = 'S8'><span style="white-space: pre"><span style="color: rgb(2, 128, 9);">%modelgeo = changeRxnBounds (modelgeo, 'EX_o2[e]', 0, 'l');</span></span></div></div><div class="inlineWrapper"><div  class = 'S8'><span style="white-space: pre"><span >modelgeo = changeRxnBounds(modelgeo, </span><span style="color: rgb(170, 4, 249);">'EX_glc_D[e]'</span><span >,-20,</span><span style="color: rgb(170, 4, 249);">'l'</span><span >);</span></span></div></div><div class="inlineWrapper"><div  class = 'S8'><span style="white-space: pre"><span >modelgeo = changeRxnBounds (modelgeo, </span><span style="color: rgb(170, 4, 249);">'EX_o2[e]'</span><span >, 0, </span><span style="color: rgb(170, 4, 249);">'l'</span><span >);</span></span></div></div><div class="inlineWrapper"><div  class = 'S8'><span style="white-space: pre"><span >modelgeo = changeObjective(modelgeo, </span><span style="color: rgb(170, 4, 249);">'DM_atp_c_'</span><span >);</span></span></div></div><div class="inlineWrapper"><div  class = 'S8'><span style="white-space: pre"><span style="color: rgb(2, 128, 9);">% WARNING: Depending on the size of the model running this function might take very long; </span></span></div></div><div class="inlineWrapper"><div  class = 'S9'><span style="white-space: pre"><span style="color: rgb(2, 128, 9);">% FBAgeo = geometricFBA (modelgeo, 'flexRel', 1e-3);</span></span></div></div></div><div  class = 'S10'><span>Display the unique fluxe</span><span>s from reactions, that are non-zero in the geometric FBA solution.</span></div><div class="CodeBlock"><div class="inlineWrapper"><div  class = 'S7'><span style="white-space: pre"><span style="color: rgb(2, 128, 9);">% for i=1:length(FBAgeo)</span></span></div></div><div class="inlineWrapper"><div  class = 'S8'><span style="white-space: pre"><span style="color: rgb(2, 128, 9);">%     if abs(FBAgeo(i)) &gt; 1e-3</span></span></div></div><div class="inlineWrapper"><div  class = 'S8'><span style="white-space: pre"><span style="color: rgb(2, 128, 9);">%         fprintf('%1.3f \t %s\n', FBAgeo(i), modelgeo.rxns{i})</span></span></div></div><div class="inlineWrapper"><div  class = 'S8'><span style="white-space: pre"><span style="color: rgb(2, 128, 9);">%     end</span></span></div></div><div class="inlineWrapper"><div  class = 'S9'><span style="white-space: pre"><span style="color: rgb(2, 128, 9);">% end</span></span></div></div></div><div  class = 'S10'><span style=' font-weight: bold;'>Troubleshooting:</span></div><div  class = 'S2'><span>When the algorithm has convergence problems, change one of the optional inputs, </span><span style=' font-family: monospace;'>flexRel</span><span>, into e.g. </span><span style=' font-family: monospace;'>1e-3</span><span>. The default is 0 when there is flexibility to flux bounds.</span></div><div  class = 'S2'><span>Enter the optional parameters as parameter name followed by parameter value, for example:  </span></div><div  class = 'S2'><span style=' font-family: monospace;'>flux = geometricFBA(model, 'epsilon', 1e-9)</span></div><h2  class = 'S6'><span style=' text-decoration: underline;'>5. Parsimonious enzyme usage Flux Balance Analysis (pFBA)</span></h2><div  class = 'S2'><span>The pFBA method was developed to achieve higher flux levels when more enzymes are required</span><span texencoding="$^5" style="vertical-align:-5px"><img src="" width="8" height="19" /></span><span>. </span></div><div  class = 'S2'><span>After performing the FBA to find the optimal value for the objective function, pFBA gets the answer of an another linear program to determine the flux distribution that minimises the total flux through all metabolic reactions in the model.</span></div><div  class = 'S2'><span style=' font-weight: bold;'>Timing:</span></div><div  class = 'S2'><span>The time to determine a pFBA solution depends on the size of the genome-scale model and is taking from </span><span mathmlencoding="&lt;math xmlns=&quot;http://www.w3.org/1998/Math/MathML&quot; display=&quot;inline&quot;&gt;&lt;mrow&gt;&lt;mo stretchy=&quot;false&quot;&gt;&amp;lt;&lt;/mo&gt;&lt;mn&gt;1&lt;/mn&gt;&lt;/mrow&gt;&lt;/math&gt;" style="vertical-align:-5px"><img src="" width="25" height="18" /></span><span> minute for a 1,000 reaction model, to 5 minutes for a model with more than 10,000 reactions.</span></div><div  class = 'S2'><span>The function is:</span></div><div class="CodeBlock"><div class="inlineWrapper"><div  class = 'S11'><span style="white-space: pre"><span style="color: rgb(2, 128, 9);">% [GeneClasses RxnClasses modelIrrevFM] = pFBA(model, varargin)</span></span></div></div></div><div  class = 'S10'><span>Where 'varagin' includes required inputs:</span></div><div class="CodeBlock"><div class="inlineWrapper"><div  class = 'S7'><span style="white-space: pre"><span style="color: rgb(2, 128, 9);">% * 'geneoption' - 0 = minimize the sum of all fluxes in the network,</span></span></div></div><div class="inlineWrapper"><div  class = 'S8'><span style="white-space: pre"><span style="color: rgb(2, 128, 9);">%                   1 = only minimize the sum of the flux through</span></span></div></div><div class="inlineWrapper"><div  class = 'S8'><span style="white-space: pre"><span style="color: rgb(2, 128, 9);">%                   gene-associated fluxes (default),</span></span></div></div><div class="inlineWrapper"><div  class = 'S8'><span style="white-space: pre"><span style="color: rgb(2, 128, 9);">%                   2 = only minimize the sum of the flux through</span></span></div></div><div class="inlineWrapper"><div  class = 'S8'><span style="white-space: pre"><span style="color: rgb(2, 128, 9);">%                   non-gene-associated fluxes</span></span></div></div><div class="inlineWrapper"><div  class = 'S8'><span style="white-space: pre"><span style="color: rgb(2, 128, 9);">%</span></span></div></div><div class="inlineWrapper"><div  class = 'S8'><span style="white-space: pre"><span style="color: rgb(2, 128, 9);">% * 'map' - map structure from readCbMap.m (no map written if empty)</span></span></div></div><div class="inlineWrapper"><div  class = 'S8'><span style="white-space: pre"><span style="color: rgb(2, 128, 9);">%</span></span></div></div><div class="inlineWrapper"><div  class = 'S8'><span style="white-space: pre"><span style="color: rgb(2, 128, 9);">% * 'mapoutname' - File Name for map</span></span></div></div><div class="inlineWrapper"><div  class = 'S8'><span style="white-space: pre"><span style="color: rgb(2, 128, 9);">%</span></span></div></div><div class="inlineWrapper"><div  class = 'S8'><span style="white-space: pre"><span style="color: rgb(2, 128, 9);">% * 'skipclass' - 0 = classify genes and reactions (default).</span></span></div></div><div class="inlineWrapper"><div  class = 'S8'><span style="white-space: pre"><span style="color: rgb(2, 128, 9);">%                 1 = Don't classify genes and reactions. Only return</span></span></div></div><div class="inlineWrapper"><div  class = 'S9'><span style="white-space: pre"><span style="color: rgb(2, 128, 9);">%                     model with the minimium flux set as upper bound.</span></span></div></div></div><div  class = 'S10'><span>Given outp</span><span>uts in this function are:</span></div><div class="CodeBlock"><div class="inlineWrapper"><div  class = 'S7'><span style="white-space: pre"><span style="color: rgb(2, 128, 9);">% OUTPUTS:</span></span></div></div><div class="inlineWrapper"><div  class = 'S8'><span style="white-space: pre"><span style="color: rgb(2, 128, 9);">% GeneClasses:  Structure with fields for each gene class</span></span></div></div><div class="inlineWrapper"><div  class = 'S8'><span style="white-space: pre"><span style="color: rgb(2, 128, 9);">% RxnsClasses:  Structure with fields for each reaction class</span></span></div></div><div class="inlineWrapper"><div  class = 'S8'><span style="white-space: pre"><span style="color: rgb(2, 128, 9);">% modelIrrevFM: Irreversible model used for minimizing flux with</span></span></div></div><div class="inlineWrapper"><div  class = 'S9'><span style="white-space: pre"><span style="color: rgb(2, 128, 9);">%               the minimum flux set as a flux upper bound</span></span></div></div></div><div  class = 'S10'><span style=' font-weight: bold;'>Calculating ATP energy production under anaerobic conditions using pFBA:</span></div><div class="CodeBlock"><div class="inlineWrapper"><div  class = 'S7'><span style="white-space: pre"><span >modelp = modelalter;</span></span></div></div><div class="inlineWrapper"><div  class = 'S8'><span style="white-space: pre"><span style="color: rgb(2, 128, 9);">%Recon3.0 model</span></span></div></div><div class="inlineWrapper"><div  class = 'S8'><span style="white-space: pre"><span >modelp = changeRxnBounds (modelp, </span><span style="color: rgb(170, 4, 249);">'EX_glc_D[e]'</span><span >, -20, </span><span style="color: rgb(170, 4, 249);">'b'</span><span >);</span></span></div></div><div class="inlineWrapper"><div  class = 'S8'><span style="white-space: pre"><span >modelp = changeRxnBounds (modelp, </span><span style="color: rgb(170, 4, 249);">'EX_o2[e]'</span><span >, 0, </span><span style="color: rgb(170, 4, 249);">'b'</span><span >);</span></span></div></div><div class="inlineWrapper"><div  class = 'S8'><span style="white-space: pre"><span >modelp = changeObjective(modelp, </span><span style="color: rgb(170, 4, 249);">'DM_atp_c_'</span><span >);</span></span></div></div><div class="inlineWrapper"><div  class = 'S8'><span style="white-space: pre"><span >[GeneClasses RxnClasses modelIrrevFM] = pFBA(modelp,</span><span style="color: rgb(14, 0, 255);">...</span></span></div></div><div class="inlineWrapper outputs"><div  class = 'S15'><span style="white-space: pre"><span >    </span><span style="color: rgb(170, 4, 249);">'geneoption'</span><span >, 0, </span><span style="color: rgb(170, 4, 249);">'skipclass'</span><span >, 1)</span></span></div><div  class = 'S13'><div class="inlineElement eoOutputWrapper embeddedOutputsTextElement" uid="AEEE734E" data-testid="output_7" data-width="420" data-height="143" data-hashorizontaloverflow="false" style="width: 450px; max-height: 261px; white-space: pre; font-style: normal; color: rgb(64, 64, 64); font-size: 12px;"><div class="textElement" style="white-space: pre; font-style: normal; color: rgb(64, 64, 64); font-size: 12px;">netFlux	fluxMeasure 	-&gt;	

GeneClasses =

     []


RxnClasses =

     []</div></div><div class="inlineElement eoOutputWrapper embeddedOutputsVariableStringElement scrollableOutput" uid="ABDDD854" data-testid="output_8" data-width="420" data-height="580" data-hashorizontaloverflow="true" style="width: 450px; max-height: 261px; white-space: pre; font-style: normal; color: rgb(64, 64, 64); font-size: 12px;"><div class="textElement" style="white-space: pre; font-style: normal; color: rgb(64, 64, 64); font-size: 12px;"><div style="white-space: pre; font-style: normal; color: rgb(64, 64, 64); font-size: 12px;"><span class="variableNameElement" style="white-space: pre; font-style: normal; color: rgb(64, 64, 64); font-size: 12px;">modelIrrevFM = <span class="headerElement" style="white-space: pre; font-style: italic; color: rgb(179, 179, 179); font-size: 12px;">struct with fields:</span></span></div><div style="white-space: pre; font-style: normal; color: rgb(64, 64, 64); font-size: 12px;">                      S: [5836×14226 double]
                   mets: {5836×1 cell}
                      b: [5836×1 double]
                 csense: [5836×1 char]
                   rxns: {14226×1 cell}
                     lb: [14226×1 double]
                     ub: [14226×1 double]
                      c: [14226×1 double]
                 osense: -1
                  genes: {2248×1 cell}
                  rules: {14226×1 cell}
             metCharges: [5836×1 int64]
            metFormulas: {5836×1 cell}
              metSmiles: {5836×1 cell}
               metNames: {5836×1 cell}
              metHMDBID: {5836×1 cell}
         metInChIString: {5836×1 cell}
              metKEGGID: {5836×1 cell}
           metPubChemID: {5836×1 cell}
            description: 'Recon3DModel.mat'
                grRules: {14226×1 cell}
             rxnGeneMat: [14226×2248 double]
    rxnConfidenceScores: [14226×1 double]
               rxnNames: {14226×1 cell}
               rxnNotes: {14226×1 cell}
           rxnECNumbers: {14226×1 cell}
          rxnReferences: {14226×1 cell}
              rxnKEGGID: {14226×1 cell}
             subSystems: {14226×1 cell}
             metCHEBIID: {5836×1 cell}
               metPdMap: {5836×1 cell}
            metReconMap: {5836×1 cell}
                modelID: 'Recon3DModel'
                 rxnCOG: {14226×1 cell}
       rxnKeggOrthology: {14226×1 cell}
            rxnReconMap: {14226×1 cell}
                version: 'Recon3D_01'
             PleaseCite: 'Brunk et al, Nat Biotech, 2018; doi:10.1038/nbt.4072'
                  match: [14226×1 double]
        reversibleModel: 0
</div></div></div></div></div></div><div  class = 'S10'><span>Display minimal fluxes of the r</span><span>eactions that are required for producing energy only from only glucose media. This is convenient with the </span><span style=' font-family: monospace;'>printFluxVector</span><span> function.</span></div><div class="CodeBlock"><div class="inlineWrapper"><div  class = 'S7'><span style="white-space: pre"><span >nonZeroFlag = 1;</span></span></div></div><div class="inlineWrapper outputs"><div  class = 'S15'><span style="white-space: pre"><span >printFluxVector(modelIrrevFM, modelIrrevFM.lb, nonZeroFlag)</span></span></div><div  class = 'S13'><div class="inlineElement eoOutputWrapper embeddedOutputsTextElement" uid="83DBC2EA" data-testid="output_9" data-width="420" data-height="45" data-hashorizontaloverflow="false" style="width: 450px; max-height: 261px; white-space: pre; font-style: normal; color: rgb(64, 64, 64); font-size: 12px;"><div class="textElement" style="white-space: pre; font-style: normal; color: rgb(64, 64, 64); font-size: 12px;">DM_atp_c_           	        1000
EX_glc_D[e]_r       	          20
netFlux             	   3.192e+04</div></div></div></div></div><h2  class = 'S6'><span style=' text-decoration: underline;'>6. Dynamic FBA</span></h2><div  class = 'S2'><span>The dynamic FBA is an extension of standard FBA that accounts for cell culture dynamics, implementing both dynamic (nonlinear programming) and static (LP) optimisation of an objective function and applying constraints to the rates of change of flux in addition to the standard FBA constraints</span><span texencoding="$^6" style="vertical-align:-5px"><img src="" width="8" height="19" /></span><span>.</span></div><div  class = 'S2'><span>The dynamic FBA metho</span><span>d implemented in this function is essentially the same as the method described by Varma A. and B. O. Palsson</span><span texencoding="$^7" style="vertical-align:-5px"><img src="" width="8" height="19" /></span><span>.</span></div><div class="CodeBlock"><div class="inlineWrapper"><div  class = 'S7'><span style="white-space: pre"><span >modeldinamic = model;</span></span></div></div><div class="inlineWrapper"><div  class = 'S8'><span style="white-space: pre"><span style="color: rgb(2, 128, 9);">% For Recon3.0 model</span></span></div></div><div class="inlineWrapper"><div  class = 'S8'><span style="white-space: pre"><span style="color: rgb(2, 128, 9);">% modeldinamic = changeRxnBounds (modeldinamic, 'EX_glc_D[e]', -20, 'l');</span></span></div></div><div class="inlineWrapper"><div  class = 'S8'><span style="white-space: pre"><span style="color: rgb(2, 128, 9);">% modeldinamic = changeRxnBounds (modeldinamic, 'EX_o2[e]', -1000, 'l');</span></span></div></div><div class="inlineWrapper"><div  class = 'S8'><span style="white-space: pre"><span style="color: rgb(2, 128, 9);">% modeldinamic = changeRxnBounds (modeldinamic, 'EX_ac[e]', -1000, 'l');</span></span></div></div><div class="inlineWrapper"><div  class = 'S8'></div></div><div class="inlineWrapper"><div  class = 'S8'><span style="white-space: pre"><span >modeldinamic = changeRxnBounds (modeldinamic, </span><span style="color: rgb(170, 4, 249);">'EX_glc_D[e]'</span><span >, -20, </span><span style="color: rgb(170, 4, 249);">'b'</span><span >);</span></span></div></div><div class="inlineWrapper"><div  class = 'S8'><span style="white-space: pre"><span >modeldinamic = changeRxnBounds (modeldinamic, </span><span style="color: rgb(170, 4, 249);">'EX_o2[e]'</span><span >, -1000, </span><span style="color: rgb(170, 4, 249);">'l'</span><span >);</span></span></div></div><div class="inlineWrapper"><div  class = 'S8'><span style="white-space: pre"><span >modeldinamic = changeRxnBounds (modeldinamic, </span><span style="color: rgb(170, 4, 249);">'EX_ac[e]'</span><span >, -1000, </span><span style="color: rgb(170, 4, 249);">'l'</span><span >);</span></span></div></div><div class="inlineWrapper"><div  class = 'S8'><span style="white-space: pre"><span style="color: rgb(2, 128, 9);">% For Recon3.0 model</span></span></div></div><div class="inlineWrapper"><div  class = 'S8'><span style="white-space: pre"><span style="color: rgb(2, 128, 9);">% smi = {'EX_glc_D[e]' 'EX_ac[e]'}; </span></span></div></div><div class="inlineWrapper"><div  class = 'S8'><span style="white-space: pre"><span >smi = {</span><span style="color: rgb(170, 4, 249);">'EX_glc_D[e]' 'EX_ac[e]'</span><span >}; </span></span></div></div><div class="inlineWrapper"><div  class = 'S8'><span style="white-space: pre"><span style="color: rgb(2, 128, 9);">% exchange reaction for substrate in environment</span></span></div></div><div class="inlineWrapper"><div  class = 'S8'></div></div><div class="inlineWrapper"><div  class = 'S8'><span style="white-space: pre"><span >smc = [10.8]; </span><span style="color: rgb(2, 128, 9);">% Glucose, Acetate concentration (all in mM)</span></span></div></div><div class="inlineWrapper"><div  class = 'S8'></div></div><div class="inlineWrapper"><div  class = 'S8'><span style="white-space: pre"><span >Xec = 0.001; </span><span style="color: rgb(2, 128, 9);">% initial biomass</span></span></div></div><div class="inlineWrapper"><div  class = 'S8'><span style="white-space: pre"><span >dt = 1.0/1000.0; </span><span style="color: rgb(2, 128, 9);">% time steps</span></span></div></div><div class="inlineWrapper"><div  class = 'S8'><span style="white-space: pre"><span >time = 1.0/dt; </span><span style="color: rgb(2, 128, 9);">% simulation time</span></span></div></div><div class="inlineWrapper"><div  class = 'S8'></div></div><div class="inlineWrapper"><div  class = 'S8'><span style="white-space: pre"><span >[concentrationMatrix, excRxnNames, timeVec,</span><span style="color: rgb(14, 0, 255);">...</span></span></div></div><div class="inlineWrapper outputs"><div  class = 'S15'><span style="white-space: pre"><span >    biomassVec] = dynamicFBA(modeldinamic, smi, smc, Xec, dt, time, smi ); </span></span></div><div  class = 'S13'><div class='variableElement' style='font-family: Menlo, Monaco, Consolas, "Courier New", monospace; font-size: 12px; '>ans = 1806</div><div class="inlineElement eoOutputWrapper embeddedOutputsTextElement" uid="90C81AC1" data-testid="output_11" data-width="420" data-height="45" data-hashorizontaloverflow="false" style="width: 450px; max-height: 261px; white-space: pre; font-style: normal; color: rgb(64, 64, 64); font-size: 12px;"><div class="textElement" style="white-space: pre; font-style: normal; color: rgb(64, 64, 64); font-size: 12px;">Step number	Biomass
Dynamic FBA analysis in progress ...
0%      [                                        ]</div></div></div></div></div><h2  class = 'S6'><span style=' text-decoration: underline;'>7. Relaxed FBA</span></h2><div  class = 'S2'><span>Find the minimal set of </span><span>relaxations on bounds and steady-state constraint to make an infeasible FBA problem feasible. </span></div><div  class = 'S2'><span>See tutorial_relaxedFBA.mlx</span></div><h2  class = 'S6'><span style=' text-decoration: underline;'>8. Flux enrichment analysis (FEA)</span></h2><div  class = 'S2'><span>The flux enrichment analysis calculates the likelihood th</span><span>at a set of fluxes would belong to a subsystem or pathway.</span></div><div  class = 'S2'><span style=' font-weight: bold;'>Timing:</span></div><div  class = 'S2'><span>The time to calcu</span><span>late the FEA is </span><span mathmlencoding="&lt;math xmlns=&quot;http://www.w3.org/1998/Math/MathML&quot; display=&quot;inline&quot;&gt;&lt;mrow&gt;&lt;mo stretchy=&quot;false&quot;&gt;&amp;lt;&lt;/mo&gt;&lt;mn&gt;1&lt;/mn&gt;&lt;/mrow&gt;&lt;/math&gt;" style="vertical-align:-5px"><img src="" width="25" height="18" /></span><span> second for any size of a model.</span></div><div class="CodeBlock"><div class="inlineWrapper"><div  class = 'S7'><span style="white-space: pre"><span >modelfea = model;</span></span></div></div><div class="inlineWrapper"><div  class = 'S8'><span style="white-space: pre"><span >res = optimizeCbModel(modelfea,</span><span style="color: rgb(170, 4, 249);">'max'</span><span >); </span></span></div></div><div class="inlineWrapper"><div  class = 'S8'><span style="white-space: pre"><span style="color: rgb(2, 128, 9);">% say you are interested in enriching the active reactions</span></span></div></div><div class="inlineWrapper"><div  class = 'S8'><span style="white-space: pre"><span >activeReactions = find(res.x) </span></span></div></div><div class="inlineWrapper"><div  class = 'S8'><span style="white-space: pre"><span style="color: rgb(2, 128, 9);">% You can also look for e.g. positive/negative/zeros flux reactions,</span></span></div></div><div class="inlineWrapper"><div  class = 'S8'><span style="white-space: pre"><span style="color: rgb(2, 128, 9);">% that depends pretty much on the question.</span></span></div></div><div class="inlineWrapper"><div  class = 'S8'><span style="white-space: pre"><span style="color: rgb(2, 128, 9);">% Now you look for the enrichement of reactions per subsystems</span></span></div></div><div class="inlineWrapper"><div  class = 'S9'><span style="white-space: pre"><span >resultCell = FEA(modelfea, activeReactions, </span><span style="color: rgb(170, 4, 249);">'subSystems'</span><span >)</span></span></div></div></div><h2  class = 'S6'><span>REFERENCES </span></h2><div  class = 'S2'><span>[1] Orth, J. D., Thiele I., and Palsson, B. Ø.  What is flux balance analysis? </span><span style=' font-style: italic;'>Nat. Biotechnol., </span><span>28(3), 245–248 (2010).</span></div><div  class = 'S2'><span>[2] Thiele, I., et al. </span><span>A community-driven global reconstruction of human metabolism. </span><span style=' font-style: italic;'>Nat. Biotechnol., </span><span>31(5), 419–425 (2013).</span></div><div  class = 'S2'><span>[3] Benyamini, T, Folger, O., Ruppin, E., Schlomi, T. Flux balance analysis accounting for metabolite dilution.</span><span style=' font-style: italic;'> Genome Biology</span><span>., 11(4):R43 (2010).</span></div><div  class = 'S2'><span>[4] Smallbone, K., and Simeonidis, E. Flux balance analysis: A geometric perspective. </span><span style=' font-style: italic;'>J Theor Biol</span><span>., 258: 311-315 (2009).</span></div><div  class = 'S2'><span>[5] Lewis, N.E., et al. Omic data from evolved E. coli are consistent with computed optimal growth from genome-scale models. </span><span style=' font-style: italic;'>Mol Syst Biol</span><span>., 6:390 (2010).</span></div><div  class = 'S2'><span>[6] Mahadevan, R., Edwards, J.S., Doyle, F.J. Dynamic Flux Balance Analysis of Diauxic Growth in Escherichia coli. </span><span style=' font-style: italic;'>Biophys J., </span><span>83(3):1331-1340 (2002).</span></div><div  class = 'S2'><span>[7] Varma A. and Palsson, B. Ø. Stoichiometric flux balance models quantitatively predict growth and metabolic by-product secretion in wild-type Escherichia coli W3110. </span><span style=' font-style: italic;'>App Environ Microbiol</span><span>., 60(10):3724-3731 (1994).</span></div><div  class = 'S2'><span>[8]. Laurent Heirendt &amp; Sylvain Arreckx, Thomas Pfau, Sebastian N. Mendoza, Anne Richelle, Almut Heinken, Hulda S. Haraldsdottir, Jacek Wachowiak, Sarah M. Keating, Vanja Vlasov, Stefania Magnusdottir, Chiam Yu Ng, German Preciat, Alise Zagare, Siu H.J. Chan, Maike K. Aurich, Catherine M. Clancy, Jennifer Modamio, John T. Sauls, Alberto Noronha, Aarash Bordbar, Benjamin Cousins, Diana C. El Assal, Luis V. Valcarcel, Inigo Apaolaza, Susan Ghaderi, Masoud Ahookhosh, Marouen Ben Guebila, Andrejs Kostromins, Nicolas Sompairac, Hoai M. Le, Ding Ma, Yuekai Sun, Lin Wang, James T. Yurkovich, Miguel A.P. Oliveira, Phan T. Vuong, Lemmer P. El Assal, Inna Kuperstein, Andrei Zinovyev, H. Scott Hinton, William A. Bryant, Francisco J. Aragon Artacho, Francisco J. Planes, Egils Stalidzans, Alejandro Maass, Santosh Vempala, Michael Hucka, Michael A. Saunders, Costas D. Maranas, Nathan E. Lewis, Thomas Sauter, Bernhard Ø. Palsson, Ines Thiele, Ronan M.T. Fleming, </span><span style=' font-weight: bold;'>Creation and analysis of biochemical constraint-based models: the COBRA Toolbox v3.0</span><span>, Nature Protocols, volume 14, pages 639–702, 2019 </span><a href = "https://doi.org/10.1038/s41596-018-0098-2"><span>doi.org/10.1038/s41596-018-0098-2</span></a><span>.</span></div><div  class = 'S2'></div>
<br>
<!-- 
##### SOURCE BEGIN #####
%% Flux Balance Analysis (FBA) and its variants
%% Author(s): *Ronan M.T. Fleming,* Leiden University;*Vanja Vlasov, LCSB, University of Luxembourg; Thomas Pfau, Systems Biology Group, LSRU, University of Luxembourg* 
%% Reviewer(s): Ines Thiele, *Catherine Clancy, National University of Ireland, Galway. Thomas Pfau, Systems Biology Group, LSRU, University of Luxembourg.* 
%% INTRODUCTION
% Flux balance analysis (FBA) evaluates the metabolic flux distribution$$^1$, 
% and is one of the most used modelling approaches for metabolic systems. 
% 
% The applications of FBA for molecular systems biology include prediction of 
% the growth rates, uptake rates, knockout lethality and product secretion. In 
% FBA, the solution space is constrained by the assumption of a steady-state, 
% under which each internal metabolite is consumed at the same rate as it is produced.
% 
% For the quantitative estimation of the metabolic fluxes, linear programming 
% (LP) can be used to solve the stoichiometric matrix for a given objective function 
% under different constraints. The constraints of the problem depict the space 
% of all eligible possibilities from which an optimal solution can be selected; 
% 
% $$\begin{array}{ll}\min\limits _{v} & \ c^{T}v\\\text{s.t.} & Sv=b,\\ & l\leq 
% v\leq u,\end{array}$$
% 
% Equation 1: Formula of standard FBA.
% 
% where $$c\in\Re^{n}$$ is a parameter vector that linearly combines one or 
% more reaction fluxes to form what is termed the objective function,  and where 
% a $$b_{i}<0$$, or  $$b_{i}>0$$, represents some fixed output, or input, of the 
% ith molecular species. $$S\in\Re^{m\times n}$$ is a stoichiometric matrix for 
% $$m$ molecular species and $$n$ reactions, and $$b$ is a vector of known metabolic 
% exchanges. The output of FBA is a particular flux distribution, $$v$ , which 
% maximises or minimises the objective function and stands between upper and lower 
% bounds, $$u$ and $$l$, respectively.
% 
% There are multiple different variants of FBA, some of which will be demonstrated 
% here:
%% 
% # *Standard FBA*
% # *Regularised FBA*
% # *Sparse FBA*
% # *Metabolite dilution FBA (mdFBA)*
% # *Geometric FBA*
% # *Parsimonious enzyme usage Flux Balance Analysis (pFBA)*
% # *Dynamic FBA*
% # *Relaxed FBA*
% # *Flux enrichment analysis (FEA)*
%% EQUIPMENT SETUP
%% Initialise The Cobra Toolbox and set the solver.
% If necessary, initialise the cobra toolbox:

 if 0
     initCobraToolbox(false) % false, as we don't want to update
 end
%% 
% For solving LP problems in a FBA analysis, certain solvers are required and 
% can be set using the |changeCobraSolver| function:

% solverOK = changeCobraSolver(solverName, solverType, printLevel, unchecked)
%% 
% The present tutorial can run with the GLPK package, which does not require 
% additional installation and configuration. Although, for the analysis of large 
% models is recommended to use the GUROBI package.
% 
% Setup the appropriate solver for the machine you are using by removing the 
% "%" (comment) sign for only the desired solver.

 changeCobraSolver('glpk','all');
% changeCobraSolver('tomlab_cplex','all');
% changeCobraASolver('ibm_cplex','all');
% changeCobraSolver ('gurobi', 'all');
%% Model Setup
% This tutorial will use the generic model of the human cellular metabolism$$ 
% ^2$, Recon 2.0. Other COBRA models, including Recon 3, may also be run with 
% this tutorial. For information on metabolites structures and reactions, and 
% to download the latest COBRA model releases, visit the Virtual Metabolic Human 
% database (VMH, <http://vmh.life). http://vmh.life).>
% 
% Before proceeding with the simulations, load the model into the workspace:

global CBTDIR
%modelFileName = 'Recon2.0model.mat';
modelFileName = 'Recon3DModel_301.mat';
modelDirectory = getDistributedModelFolder(modelFileName); %Look up the folder for the distributed Models.
modelFileName= [modelDirectory filesep modelFileName]; % Get the full path. Necessary to be sure, that the right model is loaded
if 0
    model = readCbModel(modelFileName);
else
    load('Recon3DModel_301.mat')
    model = Recon3DModel;
end
%% 
% In this tutorial we assume, that the cellular objectives include energy production 
% or optimisation of uptake rates and by-product secretion for various physiological 
% functions of the human body.
%% PROCEDURE
%% 1. Standard FBA
% Standard FBA predicts an optimal solution for a cellular objective within 
% a given set of constraints on a metabolic network (see Equation 1). Constraints 
% on the network are set by assigning limits on the uptake, consumption or production 
% of metabolites in reactions.
% 
% *Timing:*
% 
% The time to determine a FBA solution depends on the size of the genome-scale 
% model and is commonly less than a second for a medium sized model.
% 
% *Calculating maximal ATP energy production under aerobic conditions:*
% 
% For each new simulation, the original model will be copied to a new variable. 
% This preserves the constraints of the original model to perform further simulations 
% with new constraints. Additionally, this method of renaming the model avoids 
% confusion while performing multiple simulations at the same time.

modelaerobic = model;
%% 
% The ATP demand reaction, i.e., <http://vmh.life/#human/all/DM_atp_c_ |DM_atp_c_|> 
% within the model is a reaction that involves hydrolysis of ATP to ADP, Pi and 
% proton in the cytosol. 

 printRxnFormula(model, 'DM_atp_c_');
%% 
% We will set this reaction as our objective with the |'changeObjective'| command. 
% Maximising the flux through the ATP demand reaction will result in the network 
% producing a maximal amount of ATP (up to the limit of the reaction).

modelaerobic = changeObjective (modelaerobic, 'DM_atp_c_');
%% 
% The glucose and oxygen, in this case, are provided in high amounts for calculating 
% the flux through ATP demand. 
% 
% The |'changeRxnBounds'| function changes the flux constraints of the lower 
% ('|l|'), upper ('|u|'), or both the bounds ('|b|'), of the specified reaction. 
% Here, we will change the maximal uptake of glucose to 20 $\mu$mol/min/gDW and 
% of oxygen to 1000 $\mu$mol/min/gDW. The uptake of oxygen is effectively unconstrainted 
% (i.e. infinity). 

% modelaerobic = changeRxnBounds (modelaerobic, 'EX_glc_D[e]', -20, 'l'); % For Recon 3.0 uncomment these lines and 
%modelaerobic = changeRxnBounds (modelaerobic, 'EX_o2[e]', -1000, 'l'); % comment the lines below.
modelaerobic = changeRxnBounds (modelaerobic, 'EX_glc_D[e]', -20, 'l');  
modelaerobic = changeRxnBounds (modelaerobic, 'EX_o2[e]', -1000, 'l'); 
%% 
% The function |optimizeCbModel| calculates one of the optimal solutions for 
% a (maximum or minimum) objective reaction within the defined solution space. 
% In the above example, the maximal flux through the <http://vmh.life/#human/all/DM_atp_c_ 
% |DM_atp_c_|> is desired. 

FBAaerobic = optimizeCbModel (modelaerobic, 'max')
%% 
% * *Anticipated results*
%% 
% When oxygen and all carbon sources (internal and external) are provided the 
% flux through ATP demand reaction can reach its maximum rate of 1000 $\mu$mol/min/gDW. 
%% 
% * *Troubleshooting*
%% 
% If there are multiple carbon sources available in the model, it may be necessary 
% to specify more constraints in order to examine the effect of a single carbon 
% source on ATP production.  
% 
% To avoid this issue, all external carbon sources need to be closed with the 
% exception of the single carbon source of interest.

%Closing the uptake of all energy and oxygen sources
[exchBool,uptBool] = findExcRxns(model);
uptakes = model.rxns(uptBool);

% If you use Recon3.0 model, then:
% modelalter = model;
% modelalter = changeRxnBounds(modelalter, uptakes, 0, 'b');
% modelalter = changeRxnBounds(modelalter, 'EX_HC00250[e]', -1000, 'l');

% The alternative way to do that, in case you were using another large model, 
% that does not contain defined Subsystem is
% to find uptake exchange reactions with following codes:
% [selExc, selUpt] = findExcRxns(model);
% uptakes1 = model.rxns(selUpt);

% Selecting from the exchange uptake reactions those 
% which contain at least 1 carbon in the metabolites included in the reaction:
 subuptakeModel = extractSubNetwork(model, uptakes);
 hiCarbonRxns = findCarbonRxns(subuptakeModel,1);
% Closing the uptake of all the carbon sources
 modelalter = model;
 modelalter = changeRxnBounds(modelalter, hiCarbonRxns, 0, 'b');
% Closing other oxygen and energy sources. Use the following lines for recon2, or uncomment the lines below for recon3
 exoxygen = {'EX_adp[e]',   'EX_amp[e]',    'EX_atp[e]',    'EX_co2[e]',    'EX_coa[e]',    'EX_fad[e]',    'EX_fe2[e]',...
    'EX_fe3[e]',    'EX_gdp[e]',   'EX_gmp[e]',    'EX_gtp[e]',    'EX_h[e]',    'EX_h2o[e]',    'EX_h2o2[e]',...
    'EX_nad[e]',    'EX_nadp[e]',    'EX_no[e]',    'EX_no2[e]',    'EX_o2s[e]'};
modelalter = changeRxnBounds (modelalter, exoxygen, 0, 'l');
%% 
% *Calculating maximum ATP energy production under anaerobic and glucose only 
% conditions:*

modelanaerobic = modelalter;
% modelaerobic = changeRxnBounds (modelaerobic, 'EX_glc_D[e]', -20, 'l'); % For Recon 3.0 uncomment these lines and 
%modelaerobic = changeRxnBounds (modelaerobic, 'EX_o2[e]', -1000, 'l'); % comment the lines below.
modelanaerobic = changeRxnBounds(modelanaerobic, 'EX_glc_D[e]',-20,'l');
modelanaerobic = changeRxnBounds (modelanaerobic, 'EX_o2[e]', 0, 'l');
modelanaerobic = changeObjective(modelanaerobic,'DM_atp_c_');
FBAanaerob = optimizeCbModel(modelanaerobic,'max')
%% 
% * *Anticipated results*
%% 
% Compared to the aerobic condition, anaerobic condition with only glucose as 
% an energy source has reduced flux through ATP demand (82 $\mu$mol/min/gDW), 
% signifying the need to oxygen to run the oxidative phosphorylation. The results 
% are dependant on the model you are using. For Recon 3.0, under anaerobic conditions 
% with only glucose as an energy source, the flux for ATP demand is 40 $\mu$mol/min/gDW.
%% 2. Quadratically regularised FBA
% Regularised FBA calculates the optimal solution of a linear objective function, 
% and finds the unique optimal flux vector that minimises the Euclidean norm of 
% the flux vector, that is
% 
% $$\begin{array}{ll}\min\limits _{v} & \ c^{T}v + \eta ||v|| \\\text{s.t.} 
% & Sv=b,\\ & l\leq v\leq u,\\\end{array}$$
% 
% If the scalar parameter $\eta$ is small, then $c^Tv$ dominates and the same 
% optimal objective as in FBA is attained, with the advantage that a unique flux 
% vector is returned - one that minimises the the sum of the squares of the flux 
% vector.
% 
% In practice, for genome-scale, mono-scale models (e.g. not ME models) the 
% following value for the scalar parameter works fine

eta = 1e-6;
%% 
% This is a quadratic optimisation problem, so a 'QP' solver is required. The 
% COBRA Toolbox v3 comes shipped with the solver PDCO: Primal-Dual interior method 
% for Convex Objectives (further info here: <https://web.stanford.edu/group/SOL/software/pdco/ 
% link>).

[solverOK, solverInstalled] = changeCobraSolver('pdco','QP');
%% 
% Solve a quadratically regularised FBA problem:

FBAsolution = optimizeCbModel(model,'max',eta);
%% 2. Sparse FBA
% Sparse FBA calculates the optimal solution of a linear objective function 
% and finds the smallest set of reactions that can carry flux to achieve the objective. 
% Sparse FBA minimises the number of active reactions by keeping same maximal 
% objective;
% 
% $$\begin{array}{ll}\min\limits _{v} & \ ||v||_0 \\\text{s.t.} & Sv=b,\\ & 
% l\leq v\leq u,\\& c^{T}v=\rho^*\end{array}$$
% 
% Equation 2: Formula of Sparse FBA.
% 
% where the last constraint is optional and represents the requirement to satisfy 
% an optimal objective value $\rho^*$ derived from any solution to a FBA problem. 
% This approach is used to check for minimal sets of reactions that either should 
% be active or should not be active in a flux balance model that is representative 
% of a biochemical network.

% [vSparse, sparseRxnBool, essentialRxnBool]  = sparseFBA(model, osenseStr,...
%  checkMinimalSet, checkEssentialSet, zeroNormApprox)
%% 
% As an optional input, there are different appoximation types of zero-norm 
% (only available when |minNorm = 'zero'|). Default is |cappedL1|.

% Other types of zero-norm:
%  * 'cappedL1' : Capped-L1 norm
%  * 'exp'      : Exponential function
%  * 'log'      : Logarithmic function
%  * 'SCAD'     : SCAD function
%  * 'lp-'      : :math:`L_p` norm with :math:`p < 0`
%  * 'lp+'      : :math:`L_p` norm with :math:`0 < p < 1`
%  * 'l1'       : L1 norm
%  * 'all'      : try all approximations and return the best result
%% 
% *Timing:*
% 
% The time to determine a |sparseFBA()| solution depends on the size of the 
% genome-scale model and is taking from $<1$ second for a 1,000 reaction model, 
% to $<2$ seconds for a model with more than 10,000 reactions.
% 
% *Calculating maximal ATP energy production under anaerobic and glucose only 
% conditions:*

modelspar = modelalter;
% For Recon3.0 model
% modelspar = changeRxnBounds (modelspar, 'EX_glc_D[e]', -20, 'l');
% modelspar = changeRxnBounds (modelspar, 'EX_o2[e]', 0, 'l');
modelspar = changeRxnBounds(modelspar, 'EX_glc_D[e]',-20,'l');
modelspar = changeRxnBounds (modelspar, 'EX_o2[e]', 0, 'l');
modelspar = changeObjective(modelspar, 'DM_atp_c_');
[vSparse, sparseRxnBool, essentialRxnBool] = sparseFBA(modelspar, 'max');
%% 
% * *Anticipated results:*
%% 
% Commonly, a sparse FBA solution will have much smaller number of active reactions 
% compared to a standard FBA on the same model with same objective function. The 
% outputs |sparseRxnBool| and |essentialRxnBool| return vectors with 1 and 0's, 
% with sparse and essential reactions respectively.
% 
% Display the sparse flux solution, but only the non-zero fluxes. This is convenient 
% with the |printFluxVector| function.

nonZeroFlag = 1;
printFluxVector(model, vSparse, nonZeroFlag)
%% 3. Metabolite dilution flux balance analysis (mdFBA)
% This is a variant of FBA for predicting metabolic flux distributions by accounting 
% for growth-associated dilution of all metabolites in a context-dependent manner$$^3$.
% 
% A solution from the function |mdFBA| supports that all metabolites used in 
% any reaction of the solution can either be produced by the network or taken 
% up from the surrounding medium.
% 
% *Timing:*
% 
% Since this is a MIXED Integer Problem it can take a long time to solve.
% 
% *Calculating ATP energy production under aerobic condition using mdFBA:*
% 
% In this function, there is an optional output |newActives|, that represent 
% reactions that are only active in this analysis.

% The valid solution can be produced with the Recon 3.0 model
% modelmd = model;
% modelmd = changeRxnBounds(modelmd, 'EX_glc_D[e]',-20,'l');
% modelmd = changeRxnBounds (modelmd, 'EX_o2[e]', -1000, 'l');
% modelmd = changeObjective(modelmd, 'DM_atp_c_');

% [sol, newActives] = mdFBA(modelmd)
%% 
% * *Troubleshooting:*
%% 
% When a model does not have a feasible solution, add the input:  |'getInvalidSolution', 
% true|.

% clear modelmd
modelnosol = modelalter;
modelnosol = changeObjective(modelnosol, 'DM_atp_c_');
[sol, newActives] = mdFBA(modelnosol,  'getInvalidSolution', true)
%% 
% Sometimes when an FBA of a model with the same objective function and constraints 
% is run many times, or using different LP logarithm, we may get different set 
% of solutions for individual reactions. In other words, there are different sets 
% of |'FBAsolution.v'| values (fluxes of the reactions) and still get the same 
% objective function value |'f'|. That is, the opitmal solution vector is not 
% unique. This can create difficulty when investigating the changes to fluxes 
% between two different conditions. In this case a unique solution is required 
% to compare the changes to fluxes.
% 
% This issue can be solved with either of the following the methods 
%% 
% * |geometricFBA|, which provides a standard, central and reproducible solution, 
% or
% * |pFBA|, which provides a solution based on the minimal fluxes through the 
% model, and classify each gene according to how it contributes to the optimal 
% solution.
%% 4. Geometric FBA
% The geometric FBA solves the smallest frame that contains all sets of optimal 
% FBA solutions and posts a set of multiple linear programming problems$$^4$.
% 
% This FBA analysis applies iterations, where by each iteration reduces the 
% permissible solution space. After a finite number of iterations, it resolves 
% one single solution of the flux distribution.

% USAGE:
% flux = geometricFBA(model, varargin)
%% 
% *Timing:*
% 
% The time to determine a geometric FBA solution depends on the size of the 
% genome-scale model and the number of iterations. For a model with more than 
% 10,000 reactions and several iterations takes $\ge 30$minutes.
% 
% *Calculating ATP energy production under anaerobic conditions using geometric 
% FBA:*

modelgeo = modelalter;
% For Recon3.0 model
% modelgeo = changeRxnBounds (modelgeo, 'EX_glc_D[e]', -20, 'l');
%modelgeo = changeRxnBounds (modelgeo, 'EX_o2[e]', 0, 'l');
modelgeo = changeRxnBounds(modelgeo, 'EX_glc_D[e]',-20,'l');
modelgeo = changeRxnBounds (modelgeo, 'EX_o2[e]', 0, 'l');
modelgeo = changeObjective(modelgeo, 'DM_atp_c_');
% WARNING: Depending on the size of the model running this function might take very long; 
% FBAgeo = geometricFBA (modelgeo, 'flexRel', 1e-3);
%% 
% Display the unique fluxes from reactions, that are non-zero in the geometric 
% FBA solution.

% for i=1:length(FBAgeo)
%     if abs(FBAgeo(i)) > 1e-3
%         fprintf('%1.3f \t %s\n', FBAgeo(i), modelgeo.rxns{i})
%     end
% end
%% 
% *Troubleshooting:*
% 
% When the algorithm has convergence problems, change one of the optional inputs, 
% |flexRel|, into e.g. |1e-3|. The default is 0 when there is flexibility to flux 
% bounds.
% 
% Enter the optional parameters as parameter name followed by parameter value, 
% for example:  
% 
% |flux = geometricFBA(model, 'epsilon', 1e-9)|
%% 5. Parsimonious enzyme usage Flux Balance Analysis (pFBA)
% The pFBA method was developed to achieve higher flux levels when more enzymes 
% are required$$^5$. 
% 
% After performing the FBA to find the optimal value for the objective function, 
% pFBA gets the answer of an another linear program to determine the flux distribution 
% that minimises the total flux through all metabolic reactions in the model.
% 
% *Timing:*
% 
% The time to determine a pFBA solution depends on the size of the genome-scale 
% model and is taking from $<1$ minute for a 1,000 reaction model, to 5 minutes 
% for a model with more than 10,000 reactions.
% 
% The function is:

% [GeneClasses RxnClasses modelIrrevFM] = pFBA(model, varargin)
%% 
% Where 'varagin' includes required inputs:

% * 'geneoption' - 0 = minimize the sum of all fluxes in the network,
%                   1 = only minimize the sum of the flux through
%                   gene-associated fluxes (default),
%                   2 = only minimize the sum of the flux through
%                   non-gene-associated fluxes
%
% * 'map' - map structure from readCbMap.m (no map written if empty)
%
% * 'mapoutname' - File Name for map
%
% * 'skipclass' - 0 = classify genes and reactions (default).
%                 1 = Don't classify genes and reactions. Only return
%                     model with the minimium flux set as upper bound.
%% 
% Given outputs in this function are:

% OUTPUTS:
% GeneClasses:  Structure with fields for each gene class
% RxnsClasses:  Structure with fields for each reaction class
% modelIrrevFM: Irreversible model used for minimizing flux with
%               the minimum flux set as a flux upper bound
%% 
% *Calculating ATP energy production under anaerobic conditions using pFBA:*

modelp = modelalter;
%Recon3.0 model
modelp = changeRxnBounds (modelp, 'EX_glc_D[e]', -20, 'b');
modelp = changeRxnBounds (modelp, 'EX_o2[e]', 0, 'b');
modelp = changeObjective(modelp, 'DM_atp_c_');
[GeneClasses RxnClasses modelIrrevFM] = pFBA(modelp,...
    'geneoption', 0, 'skipclass', 1)
%% 
% Display minimal fluxes of the reactions that are required for producing energy 
% only from only glucose media. This is convenient with the |printFluxVector| 
% function.

nonZeroFlag = 1;
printFluxVector(modelIrrevFM, modelIrrevFM.lb, nonZeroFlag)
%% 6. Dynamic FBA
% The dynamic FBA is an extension of standard FBA that accounts for cell culture 
% dynamics, implementing both dynamic (nonlinear programming) and static (LP) 
% optimisation of an objective function and applying constraints to the rates 
% of change of flux in addition to the standard FBA constraints$$^6$.
% 
% The dynamic FBA method implemented in this function is essentially the same 
% as the method described by Varma A. and B. O. Palsson$$^7$.

modeldinamic = model;
% For Recon3.0 model
% modeldinamic = changeRxnBounds (modeldinamic, 'EX_glc_D[e]', -20, 'l');
% modeldinamic = changeRxnBounds (modeldinamic, 'EX_o2[e]', -1000, 'l');
% modeldinamic = changeRxnBounds (modeldinamic, 'EX_ac[e]', -1000, 'l');

modeldinamic = changeRxnBounds (modeldinamic, 'EX_glc_D[e]', -20, 'b');
modeldinamic = changeRxnBounds (modeldinamic, 'EX_o2[e]', -1000, 'l');
modeldinamic = changeRxnBounds (modeldinamic, 'EX_ac[e]', -1000, 'l');
% For Recon3.0 model
% smi = {'EX_glc_D[e]' 'EX_ac[e]'}; 
smi = {'EX_glc_D[e]' 'EX_ac[e]'}; 
% exchange reaction for substrate in environment

smc = [10.8]; % Glucose, Acetate concentration (all in mM)

Xec = 0.001; % initial biomass
dt = 1.0/1000.0; % time steps
time = 1.0/dt; % simulation time

[concentrationMatrix, excRxnNames, timeVec,...
    biomassVec] = dynamicFBA(modeldinamic, smi, smc, Xec, dt, time, smi ); 
%% 7. Relaxed FBA
% Find the minimal set of relaxations on bounds and steady-state constraint 
% to make an infeasible FBA problem feasible. 
% 
% See tutorial_relaxedFBA.mlx
%% 8. Flux enrichment analysis (FEA)
% The flux enrichment analysis calculates the likelihood that a set of fluxes 
% would belong to a subsystem or pathway.
% 
% *Timing:*
% 
% The time to calculate the FEA is $<1$ second for any size of a model.

modelfea = model;
res = optimizeCbModel(modelfea,'max'); 
% say you are interested in enriching the active reactions
activeReactions = find(res.x) 
% You can also look for e.g. positive/negative/zeros flux reactions,
% that depends pretty much on the question.
% Now you look for the enrichement of reactions per subsystems
resultCell = FEA(modelfea, activeReactions, 'subSystems')
%% REFERENCES 
% [1] Orth, J. D., Thiele I., and Palsson, B. Ø.  What is flux balance analysis? 
% _Nat. Biotechnol.,_ 28(3), 245–248 (2010).
% 
% [2] Thiele, I., et al. A community-driven global reconstruction of human metabolism. 
% _Nat. Biotechnol.,_ 31(5), 419–425 (2013).
% 
% [3] Benyamini, T, Folger, O., Ruppin, E., Schlomi, T. Flux balance analysis 
% accounting for metabolite dilution. _Genome Biology_., 11(4):R43 (2010).
% 
% [4] Smallbone, K., and Simeonidis, E. Flux balance analysis: A geometric perspective. 
% _J Theor Biol_., 258: 311-315 (2009).
% 
% [5] Lewis, N.E., et al. Omic data from evolved E. coli are consistent with 
% computed optimal growth from genome-scale models. _Mol Syst Biol_., 6:390 (2010).
% 
% [6] Mahadevan, R., Edwards, J.S., Doyle, F.J. Dynamic Flux Balance Analysis 
% of Diauxic Growth in Escherichia coli. _Biophys J.,_ 83(3):1331-1340 (2002).
% 
% [7] Varma A. and Palsson, B. Ø. Stoichiometric flux balance models quantitatively 
% predict growth and metabolic by-product secretion in wild-type Escherichia coli 
% W3110. _App Environ Microbiol_., 60(10):3724-3731 (1994).
% 
% [8]. Laurent Heirendt & Sylvain Arreckx, Thomas Pfau, Sebastian N. Mendoza, 
% Anne Richelle, Almut Heinken, Hulda S. Haraldsdottir, Jacek Wachowiak, Sarah 
% M. Keating, Vanja Vlasov, Stefania Magnusdottir, Chiam Yu Ng, German Preciat, 
% Alise Zagare, Siu H.J. Chan, Maike K. Aurich, Catherine M. Clancy, Jennifer 
% Modamio, John T. Sauls, Alberto Noronha, Aarash Bordbar, Benjamin Cousins, Diana 
% C. El Assal, Luis V. Valcarcel, Inigo Apaolaza, Susan Ghaderi, Masoud Ahookhosh, 
% Marouen Ben Guebila, Andrejs Kostromins, Nicolas Sompairac, Hoai M. Le, Ding 
% Ma, Yuekai Sun, Lin Wang, James T. Yurkovich, Miguel A.P. Oliveira, Phan T. 
% Vuong, Lemmer P. El Assal, Inna Kuperstein, Andrei Zinovyev, H. Scott Hinton, 
% William A. Bryant, Francisco J. Aragon Artacho, Francisco J. Planes, Egils Stalidzans, 
% Alejandro Maass, Santosh Vempala, Michael Hucka, Michael A. Saunders, Costas 
% D. Maranas, Nathan E. Lewis, Thomas Sauter, Bernhard Ø. Palsson, Ines Thiele, 
% Ronan M.T. Fleming, *Creation and analysis of biochemical constraint-based models: 
% the COBRA Toolbox v3.0*, Nature Protocols, volume 14, pages 639–702, 2019 <https://doi.org/10.1038/s41596-018-0098-2 
% doi.org/10.1038/s41596-018-0098-2>.
% 
%
##### SOURCE END #####
-->
</div></body></html>